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Abstract 


This  paper  develops  a  three  level  model  for  a  continuous  wave  diode  pumped  alkali 
laser  by  creating  rate  equations  based  on  a  three  level  system.  The  three  level  system 
consists  of  an  alkali  metal  vapor,  typically  Rb  or  Cs,  pumped  by  a  diode  from  the  2  S'  1 


state  to  the  2 Pa,  a  collisional  relaxation  from  2 Pa  to  2 Pi,  and  then  lasing  from  2 Pi 
2  2  2  2 

to  2S'i .  The  hyperfine  absorption  and  emission  cross  sections  for  these  transitions  are 
developed  in  detail.  Differential  equations  for  intra-gain  pump  attenuation  and  intra¬ 
gain  laser  growth  are  developed  in  the  fashion  done  by  Rigrod.  Using  Mathematica 
7.0,  these  differential  equations  are  solved  numerically  and  a  diode  pumped  alkali  laser 
system  is  simulated.  The  solutions  of  the  differential  equations  are  then  utilized  to 
characterize  the  inversion,  the  gain  profile,  the  output  laser  intensity,  and  the  pump 
intensity  absorption  profile  for  many  different  diode  pumped  alkali  laser  systems. 
The  results  of  the  simulation  are  compared  to  previous  experimental  results  and 
to  previous  computational  results  for  similar  systems.  The  absorption  profile  for 
the  three  level  numerical  model  is  shown  to  have  excellent  agreement  with  previous 
absorption  models.  The  lineshapes  of  the  three  level  numerical  model  are  found  to 
be  nearly  identical  to  previous  developments  excepting  those  models  assumptions. 
The  three  level  numerical  model  provides  results  closer  to  experimental  results  than 
previous  systems  and  provides  results  which  observe  effects  not  previously  modeled, 
such  as  the  effects  of  lasing  on  pump  attenuation. 
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SIMULATION  OF  A  DIODE  PUMPED  ALKALI  LASER;  A  THREE  LEVEL 

NUMERICAL  APPROACH 

I.  INTRODUCTION 

The  purpose  of  this  work  is  to  develop  a  model  for  the  propagation  of  a  Diode 
Pumped  Alkali  Laser  (DPAL)  within  a  cavity  and  to  use  computer  modeling  software 
to  implement  a  simulation  of  this  model.  The  model  will  be  developed  to  aid  in  the 
research  and  design  of  new  DPAL  systems.  A  DPAL  is  a  relatively  new  type  of  laser 
which  relies  on  laser  transitions  occurring  within  an  alkali  metal.  These  lasers  use 
electrically  driven  diodes  to  create  pump  photons,  which  are  incident  on  a  gaseous 
alkali  metal.  In  a  process  described  by  Krupke  [4],  these  photons  create  a  population 
inversion  which  leads  to  lasing.  Therefore,  DPALs  are  neither  solid  state  not  gas 
phase  laser  but  rather  a  hybrid. 

The  vast  majority  of  past  and  current  research  within  the  Air  Force  has  centered 
around  chemical  laser  systems  like  the  Chemical  Oxygen-Iodine  Laser  (COIL)  and 
Hydrogen  Fluoride  (HF)  lasers.  These  systems  offer  high  output  powers  of  at  least  a 
megawatt  or  more,  but  chemical  lasers  have  a  finite  magazine  associated  with  their 
reactants  and  usually  require  large  facilities  to  create  and  sustain  the  conditions 
needed  for  lasing.  Because  of  the  difficulty  in  deploying  chemical  laser  systems,  the 
Air  Force  is  studying  other  types  of  laser  systems  to  characterize  their  abilities  and 
their  potential  to  be  used  in  future  weapon  systems.  DPAL  systems  are  attractive 
as  one  possible  alternative  to  chemical  lasers  as  they  are  pumped  by  electrically 
driven  diodes,  and  therefore  are  capable  of  being  powered  by  a  conventional  electrical 
generator  or  an  aircraft’s  onboard  power  plant.  To  date,  DPAL  systems  are  not  well 
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characterized  compared  to  other  laser  systems.  Most  of  the  research  and  development 

on  DPAL  modeling  has  focused  on  relatively  simplistic  theoretical  models  used  to 

make  rough  estimates  and  trial  and  error  lab  characterization  of  DPAL  systems. 

In  general,  the  alkali  metal  chosen  for  a  DPAL  is  either  cesium  (Cs)  or  rubidium 

(Rb);  however,  other  alkali  metals  have  been  used  to  successfully  create  a  DPAL.  Most 

DPAL  systems  currently  are  three  level  lasers  whereby  the  ground  2Si  state  is  pumped 

to  the  excited  2Ps.  This  is  then  collisionally  relaxed  to  the  2 Pi  state  by  a  buffer  gas, 
2  2 

which  is  usually  helium  or  a  hydrocarbon  such  as  methane  or  ethane.  Photons  in 

the  2Pi  state  proceed  to  lase  to  the  ground  state.  The  2Ps  state  and  2  Pi  state  are 
2  2  2 

the  analogous  features  in  any  alkali  corresponding  to  the  well-known  doublet  in  Na. 
Other,  more  novel,  DPAL  systems  based  on  other  possible  transitions  within  the  alkali 
metals  have  been  proposed,  but  have  not  yet  been  published.  DPAL  systems  offer  a 
much  higher  stimulated  emission  cross  section  than  most  laser  species,  and  therefore, 
have  the  possibility  of  delivering  sufficient  output  intensities  and  power  needed  for  a 
weapon  system.  A  major  drawback  to  DPAL  systems  is  that  they  require  the  ground 
state  to  be  depopulated  by  at  least  half  to  create  an  inversion.  This  is  true  with 
any  three  level  laser.  To  create  the  required  inversion  demands  that  approximately 
half  of  the  atoms  in  the  2S i  be  pumped  into  one  of  the  two  upper  states.  The 
required  inversion  will  only  occur  when  fewer  atoms  are  in  the  2 Si  than  in  the  2 Pi. 
With  modern  diode  pumping  sources  this  is  achievable  albeit  somewhat  difficult. 
Currently,  the  theoretical  models  used  by  researchers  at  the  Air  Force  Institute  of 
Technology  were  developed  by  Lewis  [5]  and  Hager  [3] .  The  Lewis  and  Hager  models 
offer  many  benefits,  but  do  not  completely  characterize  a  DPAL  system.  So,  a  model 
has  been  developed  which  takes  into  account  a  greater  amount  of  parameters  and 
phenomena  than  [5]  and  [3].  Throughout  this  document,  “simulation”  refers  to  a 
computer  construct  of  a  “model”  ,  which  is  a  theoretical  construct  of  a  physical  system. 
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Simulation  and  model  are  not  interchangeable.  The  model  developed  will  provide  the 
output  power  of  a  DPAL  system  based  on  the  characteristics  of  the  pump  diodes 
and  the  lasing  cavity.  However,  this  model  will  not  be  capable  of  developing  all  the 
necessary  physics  for  the  engagement  of  a  target  with  the  DPAL  system.  Instead  this 
model  focuses  on  developing  the  output  characteristics  of  a  DPAL  system  which  then 
might  be  added  to  current  target  engagement  scenarios  or  could  be  used  as  a  research 
tool  to  direct  the  flow  of  future  DPAL  research. 
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II.  BACKGROUND 


2.1  Overview 

DPAL  systems  are  gas  electric  hybrid  lasers.  In  general,  a  gaseous  cell  of  an  alkali 
metal  has  photons  incident  on  it  from  a  narrow  banded  diode.  Alkali  atoms  are  used 
because  of  their  large  absorption  cross  sections,  well-known  properties,  and  similarity 
to  atomic  hydrogen  [5].  In  most  current  DPAL  systems,  the  diode  pumps  the  2S 1 
— y  2 Pi  transition  (the  D2  manifold).  The  incident  light  from  the  diode  must  have 
spectral  broadening  of  less  than  approximately  20  GHz  to  pump  only  the  2 Si  — >  2 Pi 
and  to  ensure  that  the  majority  of  the  pump  energy  is  absorbed  by  the  alkali.  A 
buffer  gas  is  also  present  in  the  gaseous  cell.  This  gas  serves  to  collisionally  de-excite 
atoms  from  the  state  2P a  to  the  2 Pi  state  [3].  This  transition  is  optically  forbidden, 
so  the  buffer  gas  is  required  for  a  DPAL  to  operate.  Most  often,  the  gas  chosen  is 
ethane,  methane,  helium  or  some  combination  of  the  three.  These  gases  are  generally 
selected  because  of  their  large  collisional  cross  sections  with  alkali  atoms.  The  alkali 
metals  currently  used  most  often  are  Rb  and  Cs.  These  metals  are  generally  chosen 
as  they  have  large  energy  gaps  between  the  2 Pi  and  2 Pi  states  [5].  Usually,  the 
laser  transition  occurs  between  the  2 Pi  and  the  2  Si.  The  DPAL  scheme  previously 
mentioned  is  the  most  common  method  used,  but  other  systems  have  been  theorized 
and  a  few  have  been  tested  [6].  One  such  system  uses  collisional  excitation  of  the 
alkali  to  the  2 Pa  state  by  a  collision  of  multiple  noble  gas  and  alkali  pairings.  These 
pairings  then  dissociate  after  excitation,  and  the  alkali  is  then  induced  to  lase  by  the 
same  process  aforementioned  [6].  Other  proposed  systems  include  pumping  to  upper 
states  beyond  2P ?  with  double  photon  absorption  and  creating  laser  transitions  in 
the  upper  manifold  of  the  alkali  atoms.  No  working  demonstration  of  such  a  system 
has  been  published. 
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To  date,  most  work  on  the  simulation  of  DPAL  systems  has  focused  on  modeling 
the  attenuation  of  the  diode  pump  upon  entering  the  alkali  gain  medium  and  on  de¬ 
termining  the  output  lasing  intensity  of  DPAL  systems  under  Continuous  Wave  (CW) 
conditions.  Two  of  these  simulations  were  developed  by  Hager  [3]  and  one  by  Lewis 
[5].  In  general,  the  key  to  developing  the  output  lasing  intensity  is  ascertaining  the 
number  densities  of  each  of  the  three  levels  in  the  DPAL  system,  and  then  calculating 
the  resulting  attenuation  of  the  pump  intensity.  The  intra-cavity  lasing  intensity  is 
deduced  by  the  attenuation  of  the  pump  wave  based  upon  the  conservation  of  energy. 
To  determine  these  number  densities  and  intensities,  a  simplifying  assumption  known 
as  the  quasi-two  level  system  is  used  to  some  degree  by  both  models.  However,  later 
developments  of  the  Hager  model  are  three  level  and  not  quasi-two  level. 

The  quasi-two  level  model  approximates  that  the  transition  rate  between  the  2P 3 

2 

and  2Pi  level  due  to  the  buffer  gas  is  fast  enough  that  no  other  excitation  or  de¬ 
excitation  processes  occur  to  atoms  in  the  2Ps  state.  That  is,  the  number  density 
of  the  atoms  in  the  2Ps  level  is  assumed  to  maintain  an  equilibrium  distribution  in 

relation  to  the  2 Pi  level  based  upon  statistical  mechanics.  Thus,  the  2 Pi  state  and 
2  2 

the  2  Pi  state  are  essentially  one  state.  This  approximation  is  not  completely  valid, 
but  vastly  simplifies  the  problem.  The  Lewis  model  works  under  the  quasi-two  level 
approximation,  but  attempts  to  achieve  better  fidelity  by  simulating  the  effect  of  the 
lasing  waves  in  the  cavity  after  threshold  have  been  reached.  To  actually  simulate 
the  true  effect  of  lasing,  one  must  assume  that  lasing  can  always  effect  the  number 
densities.  Without  this  effect  the  model  is  not  completely  accurate.  The  Hager  model 
assumes  that  the  number  density  along  the  alkali  cell  can  be  longitudinally  averaged, 
and  this  number  density  can  be  used  to  formulate  the  rate  equations  and  the  intra¬ 
cavity  pump  and  lasing  intensities.  This  is  known  as  Longitudinally  Averaged  Number 
Density  (LAND).  This  also,  is  not  a  wholly  accurate  approximation  and  under  certain 
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conditions  will  be  inaccurate.  Indeed,  in  most  cases  the  LAND  approximation  is  valid. 
Only  in  certain  regimes  near  threshold  does  the  LAND  approximation  begin  to  fail 
[3].  Also,  the  Hager  model  does  not  account  for  spectrally  broadband  pumping.  It  is 
assumed  that  all  pump  photons  have  the  same  frequency.  Therefore,  both  the  Lewis 
and  Hager  models  are  able  to  give  the  user  of  the  simulation  an  idea  of  how  certain 
systems  will  perform.  Both  models  offer  a  great  deal  of  insight  into  the  operation  of 
DPAL  systems  under  many  circumstances  and  regimes.  Physical  systems  cannot  be 
tested  to  sufficient  fidelity  to  replicate  the  results  of  laboratory  experiments.  Further, 
these  models  have  been  unable  to  produce  results  which  describe  a  three  level  DPAL 
system  to  the  fidelity  required  to  perform  testing  and  investigation  of  new  systems 
without  the  creation  of  an  experimental  apparatus.  Hence,  to  better  characterize 
DPAL  systems  a  more  accurate  model  of  the  DPAL  system  must  be  simulated. 

2.2  Experimental  DPAL  Development 

One  of  the  earliest  developments  of  a  laser  system  similar  to  DPAL  was  done  by 
Beach  in  2004  [1],  Beach  suggested  the  use  of  diodes  to  pump  an  optically  emitting  Rb 
or  Cs  laser  similar  to  a  diode  pumped  solid  state  laser,  but  using  a  gas  phase  alkali. 
In  this  paper,  Beach  used  a  Ti:S  laser  to  pump  a  Cs  system.  Beach  then  showed 
that  a  diode  could  also  be  used  as  a  pump  source  if  the  diode  was  able  to  emit  an 
output  sufficiently  narrow  in  frequency  to  limit  the  pumping  of  the  lasing  transition. 
To  achieve  this  a  diode’s  pump  output  must  have  a  lineshape  of  approximately  10 
GHz  at  400  Torr  of  He  to  successfully  pump  the  Cs  D2  line  [3].  In  practice,  this 
is  somewhat  difficult  to  achieve  but  diodes  have  been  successfully  used  to  pump  low 
power  systems  [10].  Beach  further  suggested  that  such  lasers  could  produce  powers 
comparable  to  solid  state  pumped  diode  lasers  in  CW  operation  [1],  A  quasi-CW 
power  of  48  W  with  a  52  percent  slope  efficiency  have  been  achieved  by  Zhdanov 
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by  pumping  with  diode  systems  [9].  It  is  believed  that  DPAL  systems  can  produce 
even  higher  output  powers  as  diode  technology  improves.  DPAL  systems  under  pulsed 
operation  have  the  potential  to  be  the  basis  of  a  megawatt  class  laser  system  [5]  So,  as 
DPAL  technology  improves  it  would  be  extremely  advantageous  and  cost  effective  to 
use  computer  modeling  to  characterize  DPAL  systems  prior  to  building  such  systems. 

2.3  Hager  Model 

The  Hager  model  can  be  either  quasi-two  level  or  three  level.  The  LAND  method 
is  utilized  to  determine  the  number  densities  used  in  the  three  level  rate  equations. 
A  two  level  model  cannot  create  the  inversion  needed  for  a  laser  [8].  The  earlier 
quasi-two  level  model  discussed  by  Hager  assumes  that  the  collisional  relaxation  rate 
between  2Ps  and  2  Pi  is  fast  enough  that  the  two  populations  are  statistically  dis¬ 
turbed  as  aforementioned.  This  model  is,  therefore,  not  a  truly  2-level  model,  but  an 
approximation  to  simplify  analysis.  The  rate  equations  are  then  developed  based  on 
these  assumptions.  It  is  assumed  that  the  laser  is  operating  in  CW  conditions.  The 
rate  equations  are  used  to  develop  the  pump  intensity  Jp  and  the  laser  intensity  II  as 
functions  of  cavity  position  via  Rigrod  analysis  [3].  The  equation  for  Ip  is  inherently 
transcendental.  The  set  of  differential  equations  is  solved  numerically. 

Hager  also  develops  the  lineshape  g]l  of  the  transitions  and  their  absorption  cross 
sections  atj .  The  gjiiy)  and  oq*  derived  include  the  effects  from  hyperfine  structure. 
The  order  of  the  subscripts  determines  whether  a  and  g  are  related  to  absorption  or 
emission.  If  i  is  first  then  it  is  absorption;  if  j  is  first  then  it  is  emission.  The  effects  of 
hyperfine  structure  on  Ip  and  Ip  are  not  examined.  As  I p  and  Ip  are  the  main  pieces 
of  information  required  to  characterize  a  DPAL  system,  Ip  and  Ip  are  the  focus  of  this 
development  [3] .  Quenching  or  collisional  de-excitation  from  levels  above  the  ground 
state  to  the  ground  state  are  not  examined.  Cavity  reflectivities  are  input  parameters 
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and  are  utilized  to  determine  loss.  The  later  three-level  model  developed  by  Hager 
assumes  that  the  number  densities  are  still  given  by  LAND,  but  that  the  2 P3  and  2  Pi 
levels  are  not  given  by  their  statistical  relation.  Each  of  these  populations  must  be 
found  independently  [3].  Both  models  developed  by  Hager  assume  single  frequency 
pumping,  which  is,  in  general,  not  valid  for  modern  diodes.  Hager  also  assumes  that 
multiple  isotopes  of  the  alkali  may  be  present  in  the  gain  medium  in  the  quasi-two 
level  and  three  level  cases.  The  Hager  model  does  simulate  the  operation  of  a  DPAL 
system  under  most  cases  well.  However,  the  Hager  model  is  not  accurate  in  all  cases 
and  can  be  improved.  Further,  in  some  systems  the  processes  not  examined  by  Hager 
can  come  to  dominate  the  operation  of  the  DPAL  system  and  must  be  handled  if  a 
complete  three  level  model  is  to  be  developed  [3]. 

2.4  Lewis  Model 

Lewis’s  model  is  a  quasi-two  level  model.  As  mentioned  previously  this  assumes 
that  the  2Ps  state  and  the  2 Pi  state  are  statistically  distributed  only.  For  simplicity, 
it  is  assumed  that  many  of  the  optical  transitions  are  much  less  likely  than  the  pump 
and  lasing  transition.  Quenching  can  be  examined,  but  is  not  because  of  a  lack  of 
data  on  the  quenching.  Several  isotopes  of  common  alkalis  are  examined  simulta¬ 
neously.  Lewis  also  handles  the  effects  of  hyperfine  structure  on  g{y)ij  and  [5]. 
Lewis  develops  the  three  level  rate  equations  and  proceeds  to  solve  them  for  the  CW 
case  based  upon  the  quasi-two  level  approximation.  The  equation  developed  for  IP 
is,  again,  transcendental.  Lewis  uses  an  analytic  solution  based  upon  the  Product 
Log  function.  This  is  the  main  benefit  of  the  Lewis  model.  Because  it  develops  an 
analytic  solution  for  I,p,  Lewis  is  able  to  provide  the  equation  for  the  pump  intensity 
as  a  function  of  position  for  any  system.  It  should  be  noted  that  the  pump  intensity 
does  not  account  for  lasing  or  the  shifting  of  population  densities  due  to  pump  in- 


tensity.  For  simplicity  it  is  assumed  that  lasing  is  not  achievable  in  general.  Lewis 
modeled  lasing  as  a  phenomena  that  would  occur  only  well  above  threshold  and  af¬ 
fected  the  pump  attenuation  globally  not  with  any  frequency  dependence  of  the  input 
photons.  Lewis  used  this  method  to  show  how  lasing  effects  the  pump  intensity  by 
drastically  increasing  the  attenuation  of  the  pump.  The  assumptions  made  to  develop 
this  model  hinder  it  from  working  accurately  in  regimes  near  threshold  and  in  cases 
when  processes  other  than  those  assumed  to  occur  come  to  dominate  the  operation 
of  the  DPAL  system  [5].  Currently,  no  model  has  handled  the  effects  of  pulsed  laser 
operation,  the  effects  of  all  possible  three  level  optical  and  collisional  transitions,  and 
the  effects  of  the  intra-cavity  pump  and  lasing  intensity  on  the  number  densities  of 
the  populations  within  the  system. 
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III.  DPAL  Theory  and  Model 


3.1  Overview 

This  chapter  outlines  the  theoretical  equations  and  chemical  kinetics  required  to 
model  a  DPAL  system.  The  three  level  system  of  alkali  atoms  used  by  DPAL  is 
explained  in  detail  and  the  hyperhne  structure  of  the  2P '3  and  2 Pi  manifolds  of  Rb 
and  Cs  is  shown.  Chemical  kinetics  are  discussed  to  develop  the  three  level  rate 

equations  and  their  dependencies.  A  lineshape  for  the  2Ps  and  2  Pi  transitions  is 

2  2 

developed  using  a  Voigt  profile.  The  lineshape  is  used  to  determine  the  stimulated 
emission  cross  sections  and  absorption  cross  sections.  The  relative  intensities  of  each 
hyperhne  transition  are  given  to  find  the  total  stimulated  emission  cross  section.  The 
rate  equations  are  then  solved  for  the  population  number  densities  known  as  N\,  N2, 
N:i.  For  a  CW  system  the  differential  equations  which  govern  attenuation  and  growth 
of  lasing  and  pump  intensities  in  the  cavity  are  given.  The  gain  coefficient  for  the 
DPAL  is  explored,  and  the  cavity  mode  spacing  is  given.  The  combination  of  these 
effects  and  equations  constitute  a  complete  picture  of  a  CW  DPAL  system. 

3.2  Three  Level  System 

The  three  level  DPAL  system  encompasses  many  possible  optical  and  collisional 

transitions.  For  this  model,  the  transitions  considered  are  given  in  Figure  1  and 

involve  the  number  densities  N:i,  N2,  and  N\  of  these  levels  respectively.  These 

transitions  are  all  of  those  between  the  2  Pa,  2  Pi,  and  2  Si  levels  which  have  population 

2  2  2 

number  densities  N3,  N2,  and  N{ .  For  the  remainder  of  this  thesis,  a  subscript  i  or 

j  denotes  the  foie  structure  energy  level,  i  and  j  can  take  on  values  of  3,  2,  or  1 

each  corresponding  to  2 Pa,  2 Pi,  and  2 Si  levels  respectively.  A31  and  A2 1  are  the 

2  2  2 

spontaneous  emission  Einstein  coefficients  (EC)  for  the  2Pa  — >  2S 1  transition  and 
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2  Pi  — »  2S'i  transition  respectively.  P13  and  P12  are  the  ECs  of  absorption  for  2 S'  1  — > 
2  2  2 

2P3  transition  and  2S'i  — >  2  Pi  transition  respectively.  P13  corresponds  to  the  pump 

transition.  P31  and  P21  are  the  ECs  for  the  stimulated  emission.  P21  corresponds 

to  the  lasing  transition.  All  /^-coefficients  are  the  rates  in  cycles  per  second  (Hz) 

of  excitation  or  de-excitation  between  an  alkali  and  some  collision  partner  (methane, 

ethane,  helium  or  an  alkali)  an  intitial  level  indicated  by  the  Erst  subscript  and  a 

final  level  indicated  by  the  second  subscript  kij  [5].  The  /^-coefficients  for  excitation 

hi 3  and  k\2  are  approximately  4  orders  of  magnitude  smaller  than  the  corresponding 

quenching  rates  k:i \  and  k2\  [5].  The  h32  and  h2 3  rates  are  the  most  important  for 

DPAL  systems  and  represent  the  collisional  relaxation  between  2Ps  and  2 Pi,  which 

allows  for  lasing  to  occur.  h3 2  and  ^23  are  approximately  4  orders  of  magnitude  larger 

than  the  quenching  rates  [5].  Transitions  to  other  atomic  states  not  listed  in  Figure 

1  can  occur,  but  they  are  much  less  likely  than  those  listed  for  a  pump  source  whose 

wavelength  closely  matches  the  wavelength  of  the  2P3  — y  2S 1  transition  for  a  given 

2  2 

alkali.  Note  that  the  transition  from  2  Pi  — y  2  Pi  is  optically  forbidden,  and  hence,  no 
ECs  are  listed  on  Figure  1.  Henceforth  “population”  and  “number  density”  are  used 
interchangeably 

From  Figure  1,  rate  equations  for  the  rate  of  change  in  time  of  each  of  the  number 
densities  iV3,  jV2,  and  Ni  can  be  constructed  in  the  same  manner  as  done  by  Lewis 
[5].  An  in  depth  treatment  is  given  in  [8].  For  the  CW  case,  these  rate  equations 
are  set  equal  to  zero  [8].  The  integrals  in  Equations  (1),  (2), and  (3)  arise  from  the 
different  possible  wavelengths  for  the  pump  and  laser  transitions  due  to  lineshape  and 
hyperhne  structure.  The  sums  over  species  are  present  to  allow  for  the  possibility  of 
different  alkali  collision  partners,  such  as  methane,  ethane,  and  helium,  being  present 
simultaneously.  Each  species  corresponds  to  a  different  collisional  excitation  and 
de-excitation  rate  and  each  must  be  handled  separately.  gp(vp)  is  the  lineshape  of 
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Figure  1.  The  three  level  chemical  kinetics  used  for  the  development  of  the  DPAL 
model.  N i  corresponds  to  the  population  in  the  2Si  state.  N 2  corresponds  to  the 
population  in  the  2  Pi  state.  N 3  corresponds  to  the  population  in  the  2Pa  state.  The 
pump  transition  corresponds  to  P13.  The  laser  transition  corresponds  to  P21  • 
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the  input  pump  as  a  function  of  up.  The  concentration  of  a  given  collision  partner  is 
M Speaes .  Also,  the  excitation  rate  coeffcients  k\3  and  k\ 2  are  neglected  because  of  their 
extremely  small  size  under  normal  temperatures  for  operation  of  a  DPAL  system  (less 
than  approximately  1000  K)  and  for  most  collision  partners  typically  used  (methane, 
ethane,  and  helium).  However,  if  one  were  to  use  another  gas  as  a  collision  partner 
these  terms  might  need  to  be  reinserted.  It  is  suggested  that  the  reader  use  the  list 
of  symbols  provided  at  the  outset  of  this  document  while  reading  these  equations  to 
avoid  confusion,  h  is  Planck’s  constant. 
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Equations  (1),  (2), and  (3)  may  then  be  set  equal  to  zero  in  the  CW  case.  The 
CW  case  also  eliminates  the  dependence  on  time  in  Ip(u,z,t)  and  Ii(u,z,t).  These 
equations  constitute  two  independent  equations,  not  three,  as  all  of  the  rates  between 
the  three  levels  are  present  in  two  of  the  equations  simultaneously.  Hence,  these 
equations  do  not  form  an  independent  set  of  three  equations,  but  only  two  dependent 
equations.  Therefore,  to  solve  for  the  corresponding  population  number  densities  Ni, 
N2,  and  N3,  which  is  the  ultimate  goal  of  this  development,  another  equation  must  be 
added.  This  is  Equation  (4),  where  Nt  represents  the  total  number  density  of  alkali 
in  the  system. 


Nt  =  Ni  +  N2  +  N3 


(4) 


3.3  Cross  sections  and  Lineshapes  for  DPALs 

With  three  independent  equations,  Ni,  N2,  and  N3  may  be  solved  unambiguously. 
This  is  performed  via  Mathematica  in  Appendix  A.  With  this  level  of  detial  in  the 
rate  equations,  this  is  a  non-trivial  task  without  a  computational  aide.  Armed  with 
N\ ,  Ni,  and  N3  the  individual  terms  within  the  rate  equations  must  now  be  developed 
to  actually  compute  the  populations  A1;  N2,  and  N3.  The  stimulated  emission  cross 
section  aJt  for  a  DPAL  system  is  given  by  Equations  (5)  and  (6)  where  cr^-  is  the 
corresponding  absorption  cross  section.  If  j  and  %  correspond  to  N3  and  N{ ,  the 
pump  frequency  up  is  used,  and  if  j  and  i  correspond  to  N2  and  W,  the  pump 
frequency  zy  is  used.  The  lineshape  gjiiy ),  the  natural  abundance  of  hyperhne  states 
fji(F",F'),  and  relative  intensity  of  the  hyperhne  transitions  S(F" ,  F' ,iso)  used  all 
follow  the  same  employment  scheme  for  i  and  j  as  atj .  It  is  of  note,  that  while  the 
lasing  frequency  is  a  single  value  in  CW  single  mode  operation  as  per  [8],  the  pump 
frequency  up  has  a  lineshape  associated  with  the  output  of  the  diode.  This  requires 
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that  while  cr2i  (/y)  has  a  single  value,  <t31  (up)  is  a  continuous  function.  This  point  is 
easily  overlooked  and  is  important  to  accurately  model  a  broadband  pumped  DPAL 
system.  The  functional  dependence  on  a  particluar  isotope  is  shown  as  {iso)  and 
relates  to  the  particular  alkali  isotope  being  used.  In  Equation  (5),  iso  denotes  the 
particular  isotope  species  and  its  particular  gji(v),  S(F" ,  F' ,iso),  and  fji{F",iso ). 
fiSO  is  the  natural  abundance  of  a  particular  isotope  as  a  fraction  out  of  one. 

2 

<Jij(l')  =  Y  Y  Y  —^^AjiS{iso^F'^F')9ji{v)fji{F'\iso)fiso  (5) 

iso  F”{iso)  F'(iso) 


GjiW)  =  —Gij{v)  (6) 

9j 

The  lineshape  gjj  is  a  convolution  of  the  Doppler  broadening  (a  Gaussian),  the 
pressure  broadening  (a  Lorentzian),  and  the  lifetime  broadening  (a  Lorentzian)  [8] . 
The  lifetime  and  pressure  broadening  are  known  as  the  homogenous  broadening  and 
are  given  by  A vhji(iso)  in  Equation  7.  The  Doppler  broadening  is  A vdji(iso)  given 
by  Equation  (8)  [8].  Using  these  quantities,  the  Voigt  lineshape  gji(v)  is  given  by 
Equation  (9)  from  [8].  kb  is  Boltzmann’s  constant. 
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A udji(F" ,  Fr,  iso)  depends  on  the  particular  hyperhne  line  to  which  A udji(F" ,  F iso) 
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corresponds  as  each  hyperfine  transition  has  its  own  associate  FWHM.  The  FWHM 
depends  upon  F",  the  original  transition  state  in  either  N\ ,  jV2,  or  N3,  and  F',  the 
final  state  in  Afi,  N2,  or  N3  [5].  Note,  vhyji(F' ,  F" ,iso)  is  the  hyperfine  frequency, 
in  Hz ,  for  a  transition  from  one  state  to  another  for  a  particular  F',  F" ,  and  iso¬ 
tope.  vhyji(F' ,  F" ,iso)  is  not  the  splitting  between  hyperfine  levels  within  a  single 
fine  level.  That  is  ^hyspHtj(F” ,  iso).  Figures  2,  3,  4,  5,  6,  and  7  show  the  transitions 
between  each  of  the  hyperfine  states  of  the  D1  and  D2  transitions  for  Cs,  85Rb,  and 
87 Rb.  The  selection  rule  for  F"  to  F'  is  A F—  +1,-1,  or  0  with  F"  =  0  to  F'  =  0 
being  forbidden.  Figures  2,  3,  4,  5,  6,  and  7  are  taken  from  [5]. 

In  practice,  Equation  (9)  is  difficult  to  use.  An  approximation  using  an  error 
function  can  be  made  as  is  employed  as  per  [7] .  The  intermediate  quantities  given  by 
Equations  (10),  (11),  and  (12)  are  developed  in  this  process.  In  Equation  (12),  the  i 
in  front  of  Uji(v,  F" ,  F' ,iso)  is  the  square  root  of  negative  one.  For  the  remainder  of 
this  thesis,  Equation  (13)  will  be  used  for  gji(v,  F" ,  F' ,iso),  not  Equation  (9). 


dji(v,  F",  F',  iso)  =  Log( 2)  2 


A  vhji  (iso) 
A  vdji(iso) 


(10) 


Uji(u,  F" ,  F' ,  iso) 


2Log(2) 


1  iy_  -  v hyji(F",F',iso )) 
A odji(F",  F' ,  iso) 


(11) 


Zji(v,  F" ,  F',  iso) 


=  a,ji(v ,  F" ,  F' ,  iso)  +  %  =t=  Uji(v ,  F",  F' ,  iso) 


(12) 


gji{ui  F"  1 F1,  iso)  = 


4Lo5(2)1/2 


7 r 


Aisdji(F ",  F\  iso 


Re(eZji(l',F"’F'’is°)2 ErfC(Zji(v,  F" ,  F' ,  iso))) 


(13) 

To  calculate  the  lineshape,  the  relative  strength  S (species,  F" ,  F' ,iso)  of  a  tran¬ 
sition  between  hyperfine  levels  must  be  known  [5].  The  relative  intensities  of  each  of 
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Figure  2.  The  D1  manifold  for  133Cs  without  possible  transitions  listed.  The  frequency 
spacing  of  each  hyperflne  state  is  listed  [5]. 
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Figure  3.  The  D2  manifold  for  133Cs  without  possible  transitions  listed.  The  frequency 
spacing  of  each  hyperflne  state  is  listed  [5]. 
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Figure  4.  The  D1  manifold  for  85Rb  without  possible  transitions  listed.  The  frequency 
spacing  of  each  hyperflne  state  is  listed  [5]. 
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Figure  5.  The  D2  manifold  for  85Rb  without  possible  transitions  listed.  The  frequency 
spacing  of  each  hyperflne  state  is  listed  [5]. 
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Figure  6.  The  D1  manifold  for  s'Rb  without  possible  transitions  listed.  The  frequency 
spacing  of  each  hyperflne  state  is  listed  [5]. 
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Figure  7.  The  D2  manifold  for  s'Rb  without  possible  transitions  listed.  The  frequency 
spacing  of  each  hyperfine  state  is  listed  [5]. 
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the  transitions  for  the  D1  and  D2  transitions  may  be  found  in  Tables  1,  2,  and  3. 
The  data  in  Tables  1,  2,  and  3  comes  from  [5]. 

The  relative  natural  abundance  of  each  hyperhne  state  as  a  fraction  of  one  is 
fji(F",iso).  fji(F",iso)  is  given  by  Equation  (14).  The  relative  natural  abundance 
is  the  fraction  of  atoms  in  each  of  the  initial  hyperhne  states  F" .  It  is  of  note  that 
Vhy split j(F",  iso)  is  not  equal  to  vhyji(F" ,  F' ,  iso) .  Vhyspiit^F" ,iso)  is  the  frequency  in 
Hz  of  the  splitting  between  a  particular  hyperhne  level  and  the  center  frequency  of 
the  fine  structure  level  which  the  hyperhne  level  is  in.  While,  v,hyji(F",F,:iso )  is  the 
frequency  in  Hz  of  a  F"  — >  F1  hyperhne  transition  between  two  fine  structure  levels. 
Simply  put,  Vhy  split  3(F"  ,iso)  is  an  energy  splitting  in  Hz  given  on  Figures  2,  3,  4,  5, 
6,  and  7  and  vhy..(F" ,  F' ,iso)  is  the  absolute  difference  between  two  hyperhne  levels 
in  different  fine  structure  levels. 


fji(F",  iso) 


(2  F"  +  1  )exp( 


h(Vji  Vhysplitz  (.F 

Wr 


+  l)exp( 

iso  F" 


split  z  (F" , iso) 

Wr 


3.4  Alkali  and  Collision  Partner  Properties 


(14) 


Several  other  important  material  quantities  are  needed  to  calculate  Equation  (6). 
Table  4  provides  these  quantities  for  Cs  and  Rb. 

From  the  numbers  in  Table  4  and  using  equations  from  [8],  several  more  needed 
quantities  are  developed.  These  include  Bji  and  kji. 


Bji(v,iso) 


c3Aji(iso ) 
n387r/iu3 


(15) 


Bij(v,  iso) 


— Bji(v,iso ) 
9i 


(16) 
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Table  1.  S (species,  F” ,  F' ,  iso)  for  133Cs  for  Hyperfine  Structure 


D1  Manifold 

D2  Manifold 

F"  -)•  F'  or  2Si  -)■  2P i 

2  2 

S (species,  F" ,  F',  iso ) 

F"  -)•  F'  or  2Pi  -)■  2Ps 

2  2 

S(species,  F",  F',  iso ) 

3  — »  3 

1/4 

3  — >■  2 

5/14 

3—^4 

3/4 

3  ->  3 

3/8 

4  — )■  3 

7/12 

3  — )■  4 

15/56 

4  — >■  4 

5/12 

4  — )■  3 

7/72 

- 

- 

4  — >■  4 

7/24 

- 

- 

4  — )■  5 

7/18 

Table  2.  S (species,  F" ,  F' ,  iso)  for  85Rb  for  Hyperfine  Structure 


D1  Manifold 

D2  Manifold 

F"  -)■  F'  or  2Si  -)■  2Fi 

2  2 

S (species,  F" ,  F',  iso ) 

F"  -)•  F'  or  2Pi  ->■  2Fs 

2  2 

S(species,  F" ,  F’ ,  iso ) 

2  — >■  2 

2/9 

2  ->■  1 

3/10 

2  — )■  3 

7/9 

2  — >■  2 

7/18 

3  — )■  2 

5/9 

2  — )■  3 

14/45 

3—^3 

4/9 

3  — >■  2 

5/63 

- 

- 

3—^3 

5/18 

- 

- 

3  — )■  4 

9/14 

Table  3.  S (species,  F" ,  F' ,  iso)  for  8'Rb  for  Hyperfine  Structure 


D1  Manifold 

D2  Manifold 

F"  ->•  F’  or  2Si  -)•  2Fi 

2  2 

S (species,  F",  F' ,  iso ) 

F"  -f  F’  or  2 Pi  ->•  2Fa 

2  2 

S(species,  F" ,  F' ,  iso ) 

1  ->■  1 

1/6 

1  0 

1/6 

1  -)•  2 

5/6 

1  ->■  1 

5/12 

2  ->•  1 

1/2 

1  2 

5/12 

2  — )■  2 

1/2 

2  ->•  1 

1/20 

- 

- 

2  — >■  2 

1/4 

- 

- 

2  — )■  3 

7/10 
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Table  4.  133Cs,  85Rb,  and  87Rb  Material  Properties 


133  Cs 

85  Rb 

87Rb 

fiso 

1 

0.2783  [3] 

0.7217[3] 

i 

7/2(5] 

5/2(3] 

3/2(3] 

A2i(Hz) 

n/a 

27.679X10-9 

27.7x10-9  [5] 

A31(Hz) 

n/a 

26.2348x10-9 

26.24x10-9  [5] 

93 

4  [5] 

4  [5] 

4  [5] 

92 

2  (5] 

2(5] 

2(5] 

9i 

2  (5] 

2(5] 

2(5] 

7 He  for  Dl(MH  z/Torr ) 

26.21  [5] 

18.9  ±0.2  [9] 

18.9  ±0.2  [9] 

7 He  f°r  D2(MH  z/Torr ) 

23.50  [5] 

20.0  ±0.4  [9] 

20.0  ±0.4  [9] 

5 He  for  Dl(MH z/Torr) 

4.46  [5] 

4.71  ±0.04  [9] 

4.71  ±0.04  [9] 

5 He  f°r  D2(MH z/Torr) 

0.75  [5] 

0.37  ±0.06  [9] 

0.37  ±0.06  [9] 

1cha  for  Dl(MH z/Torr) 

29.1  [5] 

29.1  ±0.8  [9] 

29.1  ±0.8  [9] 

7 Ch4  for  D2(MH z/Torr) 

26.2  [5] 

26.2  ±0.6  [9] 

26.2  ±0.6  [9] 

5cha  for  D\[MH z/Torr) 

-7.92  [5] 

-7.9  ±0.1  [9] 

-7.9  ±0.1  [9] 

bcnA  for  D2{MH z/Torr) 

-7.0  [5] 

-7.0  ±0.1  [9] 

-7.0  ±0.2  [9] 

k32  (cm3s~1) 

? 

3.16  x  10~10  [2] 

3.16  x  10“10  [2] 

k3 1 

o  [2] 

0  [2] 

0  [2] 

k-2 1 

0  [2] 

0(2] 

0  [2] 

kij{iso)  =  kji{iso ) 


9i  / 

— exp( 

9j 


hvji 

kbT 


3.5  Lasing  and  Pump  Intensities 


(17) 


To  find  the  output  lasing  intensity  one  must  first  determine  the  intra-cavity  lasing 
intensity  and  the  intra-cavity  pump  intensity  as  functions  of  up  and  time.  A  nearly 
accurate  treatment  of  this  is  done  by  [8]  and  is  known  as  Rigrod  analysis.  Rigrod 
analysis  involves  the  propagation  in  the  cavity  of  a  pump  wave,  Ip,  a.  forward  trav¬ 
eling  lasing  wave,  I+,  and  a  backward  traveling  lasing  wave,  /_  [8].  The  saturation 
intensity  Isat  is  also  required  to  perform  Rigrod  analysis  and  is  given  by  Equation 
(18).  Equations  (18),  (19),  (20),  and  (21)  are  the  initial  differential  equations  for  the 
attenuation  or  growth  of  intra-cavity  pump  and  lasing  waves.  Notice  that  Equations 
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(19),  (20),  and  (21)  are  transcendental.  Also,  remember  that  Ah,  N2,  and  N3  depend 
on  Ip,  I+,  and  This  dependence  is  extremely  lengthy  and  is  shown  in  Appendix 
A.  These  two  facts  greatly  increase  the  solution  difficulty  of  these  differential  equa¬ 
tions.  Equations  (19),  (20),  and  (21)  require  numerical  differential  equation  solution 
techniques  to  solve.  Equations  (20)  and  (21)  appear  without  the  terms  which  go  as 
the  inverse  of  I+  plus  /_  divided  by  Isat  because  the  population  densities  (Ah,  AT2, 
and  N3)  maintain  their  dependence  on  Ip,  I+,  and  /_  [8]. 


Isatiyp)  — 


hun 


E 


1 


031  (z'p)  .  Trad(isO ) 


fiso ) 


(18) 


dlp{z,  i>i,  Up) 
dz 


g  3 

031  iyp)  [H3(z,  Vu  Vp) - Ni(z,  Vh  Vp)\Ip(z,  Vh  Vp) 

9i 


(19) 


dl+(z,ui,up ) 
dz 


92 

021  (n;) [Ah {z,  Vi,  Up) - Nx(z,  uh  up)\I+{z,  uu  up) 

9 1 


(20) 


(z,ui,up)  =  -a2i{ui)[N2(z,ul,up)  -  —Ni(z,ui,vp)\I+(z,vi,up)  (21) 

dz  9i 

The  crux  of  the  simulation  lies  in  solving  Equations  (19),  (20),  and  (21).  The 
solution  to  these  differential  equations  requires  three  boundary  conditions.  Initial 
conditions  are  not  used  as  the  solution  is  time-independent  for  the  CW  case.  Equa¬ 
tions  (22),  (23),  and  (24)  provide  these  boundary  conditions  as  given  by  [8]  for  Rigrod 
analysis,  where  ud  is  the  line  center  of  the  diode’s  spectral  profile  in  Hz.  Zq  is  the  po¬ 
sition  of  the  beginning  of  the  gain  cell  and  zj  is  the  position  of  the  end  of  the  gain  cell. 
Note,  that  z  denotes  the  functional  dependence  of  an  equation  upon  longitundinal 
position  in  m  in  the  alkali  gain  cell. 

Ip(Zo,ul,up)  =  Ip0^p\  (22) 

93l{Vd) 
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i+Oo,  Vi,  Vp)  =  RiTgI^(z0,  Vi,  Up) 


(23) 


I-(zf,  uh  Up)  =  R2Tp+(zf,  uh  Up)  (24) 

Other  important  and  oft  used  benchmarks  for  the  system  are  the  gain  coefficient 
7„  given  by  Equation  (25),  the  loss  coefficient  a  given  by  Equation  (26),  and  the 
cavity  mode  spacing  A Ufsr  given  by  Equation  (27).  Under  CW  operation,  the  gain 
coefficient  will  equal  the  loss  coefficient  once  threshold  is  reached  [8].  This  forces  the 
gain  to  operate  at  the  loss.  The  cavity  mode  spacing  and  the  gain  coefficient  are 
important  as  they  determine  what  the  operating  frequency  of  the  laser  {yi)  will  be 
as  the  mode  with  the  greatest  gain  coefficient  will  be  selected  preferentially  over  all 
other  modes  [8].  By  the  time  threshold  is  reached,  all  of  the  power  in  the  laser  will  be 
focused  in  this  single  mode.  This  is  the  idea  of  single  mode  operation.  This  is  only 
true  for  lasers  operating  in  a  single  TEM  mode  and  not  multiple  modes  [8].  n  is  the 
index  of  refraction,  c  is  the  speed  of  light. 


7 v{z,uhup)  =  r/i+(~o \,  vi,Vp)  — 1 -  =  0-21(i>j)(jV2(z,  ut,  Up)  -  —Nl(z,uuup))  (25) 

dz  I+{z0,Ui,Up )  g  l 


a=wgln{Rk] 


(26) 


Vq+ 1  ~Vq  = 


'2jl  (l  m  /  f  ror 

These  equations  constitute  a  completely-modeled  DPAL  system. 


(27) 
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IV.  Simulation  Description 


4.1  Overview 

The  simulation  of  the  model  developed  in  Chapter  III  is  performed  in  Mathematica 
7.0  for  Windows  XP  on  a  2.0  GHz  AMD  processor.  The  simulation  reads  in  a  list  of 
input  parameters,  then  uses  all  of  the  equations  developed  in  Chapter  III  to  simulate 
a  DPAL  system.  Then  using  the  differential  equation  solving,  data  analysis,  and  visu¬ 
alizations  packages  in  Mathematica  the  system  is  characterized  as  detailed  in  Chap¬ 
ter  III.  The  outputs  required  of  the  system  are  Ip(z,  zy,  up),  I+(z,  zy,  up),  J_(z,  ui,  up), 
7 u{z,  zy,  up),  a,  Ni(z,  uh  up),  N2(z,  uh  up),  N3(z,  ut,  up ),  and  the  output  laser  intensity, 
lout ,  which  is  the  output  at  the  single  laser  frequency,  zy 

4.2  Assumptions 

Without  simplifying  assumptions  simulating  the  DPAL  model  discussed  in  Chap¬ 
ter  III  would  not  be  practical  as  a  single  threaded  process  on  a  desktop  computer. 
Assumptions  were  made  to  simplify  the  problem,  but  these  assumptions  were  chosen 
such  that  a  high  degree  of  fidelity  is  maintained  in  the  simulation’s  outputs.  The 
DPAL  system  is  assumed  to  be  CW  in  its  pump  and  response.  Hence,  all  time  depen¬ 
dence  is  eliminated.  For  any  system  which  runs  for  longer  than  approximately  one 
ms  this  is  an  adequate  assumption  because  the  population  densities  will  reach  their 
equilibrium  values  within  that  time.  Only  transitions  within  the  three  level  system 
shown  in  Figure  1  are  assumed  to  occur.  All  other  transitions  within  the  alkali  atom 
are  ignored.  This  assumption  is  somewhat  valid  as  the  absorption  and  emission  cross 
section  for  the  2Si  — >  2P3  and  2  Pi  — >  2 S' i  are  much  larger  (approximately  104  times 
larger)  than  any  of  the  other  transitions. 

It  is  assumed  that  no  forbidden  hyperhne  transitions  occur  within  or  between  the 


sublevels.  The  effects  of  the  hyperfine  structure  beyond  its  effects  on  the  lineshape 
gjiiy )  are  neglected.  The  hyperfine  structure  in  a  physical  system  does  play  a  part  in 
the  rate  between  each  sublevel.  Hence,  a  fully  complete  model  would  need  to  create 
rate  equations,  such  as  Equations  (1) , (2) ,  and  (3),  for  each  possible  (including  the 
forbidden  transitions)  hyperfine  transition.  The  populations  of  each  hyperfine  level 
would  also  have  to  be  independent  of  one  another.  Transitions  between  hyperfine 
levels  due  to  collisions,  absorption,  and  emission  would  also  have  to  be  tracked.  Fur¬ 
ther,  the  lasing  intensity  would  not  be  single  quantity  but  a  set  of  intensities  each 
associated  with  a  transition  on  the  D1  manifold.  Each  intensity  would  require  an  in¬ 
dividual  plus  wave  and  minus  wave  differential  equation  like  Equations  (20)  and  (21). 
Thus,  the  system  would  contain  not  three  coupled  non-linear  ODEs,  but  instead  13. 
Also,  the  collisional  excitation  and  de-excitation  rates  between  hyperfine  levels  are 
unknown.  Adding  this  complexity  to  the  model  is  both  extremely  complicated  and 
computationally  difficult  without  the  use  of  multi-threading.  The  laser  modeled  is  as¬ 
sumed  to  operate  only  in  the  TEM(0,0)  mode  and  to  not  operate  in  any  other  modes. 
The  mode  volume  is  assumed  to  be  completely  filled.  The  energy  lost  because  of 
unfilled  mode  volume  is  not  considered.  This  assumption  is  only  somewhat  valid,  but 
it  is  difficult  to  simulate  the  effects  of  partially  filling  the  mode  volume  of  a  DPAL 
because  these  effects  are  not  yet  fully  characterized  in  the  literature. 

4.3  Simulation  Input  Parameters 

The  inputs  to  the  simulation  are  listed  in  Table  5.  Each  parameter  is  given  with 
its  normal  range  of  values.  The  input  parameters  are  not  physical  constants,  but 
are  rather  variables  which  are  dynamic  between  different  DPAL  systems.  Collisional 
relaxation  rate  coefficients  such  as  are  not  considered  input  parameters.  These 
are  considered  physical  constants  of  the  system  and  are  not  input  by  the  user.  The 
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simulation  has  the  ability  to  include  more  collisional  relaxation  rates  than  those  for 
methane  and  helium,  but  does  not  currently.  Appendix  B  shows  a  list  of  the  physical 
parameters  used  by  the  code.  It  is  of  note  that  the  three  level  numerical  model  takes 
all  inputs  in  meters,  kilograms,  and  seconds  (mks)  and  the  derived  units  of  the  mks 
system.  Units  with  metric  prefixes  other  than  kilograms  should  not  be  used.  Table  5 
presents  the  inputs  with  prefixes  for  simplicity. 


Table  5.  Simulation  Input  Parameters 


Input  Parameter 

LInits  required 

Normal  Range 

Symbol 

Temperature  of  Cell 

K 

250-2000 

T 

Temperature  of  Methane 

K 

250-2000 

Tmeth 

Temperature  of  Helium 

K 

250-2000 

THe 

Temperature  of  Alkali 

K 

250-2000 

Talk 

Partial  Pressure  of  Methane 

kPa 

20-2000 

Mmeth 

Partial  Pressure  of  Helium 

kPa 

20-2000 

MHe 

Partial  Pressure  of  Alkali 

kPa 

2-200 

Malk 

Total  Alkali  Number  Density 

m 

1018-1024 

Nt 

Length  of  gain  medium 

m 

0.005-0.5 

h 

Transmission  of  windows  around  gain 

LI  nit  less 

0.95-1.0 

Tg 

High  Reflector  Reflectivity 

LI  nit  less 

0.95-1.0 

Ri 

Output  Coupler  Reflectivity 

LI  nit  less 

0. 2-1.0 

R-2 

Initial  Pump  Intensity  at  linecenter 

Wcm~'2 

102  -  106 

IpO 

Diode  Pump  Line  Center  frequency 

THz 

384.0  -  384.4 

Vd 

Diode  Pump  FWHM  frequency 

GHz 

o 

1 

1 

o 

to 

^ Pfwhm 

Space  Between  Cavity  Mirrors 

m 

0.01-1.0 

^mirror 

4.4  Simulation  Outline 

The  simulation  begins  by  reading  the  input  parameters  in  Table  5.  Then  the  pro¬ 
gram  develops  each  equation  listed  in  Chapter  III  from  Equation  (1)  to  (27)  except 
Equation  (9),  to  which  an  approximation  was  used  as  detailed  in  Chapter  III.  These 
equations  are  supplemented  by  the  physical  parameters  of  the  system  given  in  Tables 
1-4  and  Figures  2-7.  During  this  development  the  assumptions  listed  above  are 
made  by  the  simulation  or  were  made  during  the  simulation’s  design.  The  code  then 
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constitutes  a  fully-developed  model  of  a  DPAL  system.  If  the  differntial  equations 
(19),  (20),  and  (21)  can  be  solved  numerically  for  this  system,  then  all  of  the  pa¬ 
rameters  needed  to  characterize  the  system  can  be  developed  from  that  solution.  In 
practice,  the  solution  of  Equations  (19),  (20),  and  (21)  is  difficult  and  requires  the 
use  of  Mathematical s  innate  numerical  differential  equation  solver.  The  differential 
equation  solver  is  utilized  with  a  shooting  method  and  the  boundary  conditions  of 
Equations  (22),  (23),  and  (24).  To  solve  Equations  (19),  (20),  and  (21)  requires  that 
a  shooting  method  be  performed  for  different  initial  values  for  I+  and  /_  at  zq.  Each 
shot  is  then  solved  by  the  numerical  differential  equation  solver  for  a  solution  to  the 
differential  equations  in  z  or  position.  An  algorithm  then  assigns  the  best  starting  I+ 
and  /_  values  based  upon  the  occurrence  of  gain  above  loss  and  on  the  degree  of  ad¬ 
herence  to  Equation  (24).  This  solution  architecture  is  then  utilized  for  a  select  set  of 
pump  frequencies  spanning  three  times  the  pump’s  full  width  at  half  max  (FWHM). 
The  solutions  which  best  match  the  boundary  condition  for  each  discrete  pump  fre¬ 
quency  are  then  interpolated  between  by  Mathematica' s  data  analysis  software  and 
are  plotted.  This  architecture  is  detailed  in  Figure  8.  The  complete  Mathematica 
notebook  can  be  found  in  Appendix  B. 

4.5  Simulation  Outputs 

The  outputs  of  the  simulation  are  the  way  in  which  the  model  is  both  com¬ 
pared  to  other  models  such  as  [2], [3],  and  [5],  and  then,  the  outputs  are  compared 
to  experimental  results.  To  achieve  these  comparisons  several  different  outputs  are 
provided.  A  profile  of  the  pump  intensity  Ip(z,up)  is  shown  throughout  the  entire 
cell.  I+(z,  up,  zy)  and  I-(z,  up,  zy)  are  plotted.  The  population  densities,  N\  (z,  up), 
N2(z,  Up),  and  N3(z,  up )  are  plotted.  The  output  lasing  intensity  (lout)  is  determined. 
The  simulation  selects  the  lasing  mode  with  the  greatest  gain  and  outputs  its  fre- 
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Physical  Parameters  and  Constants 


Lineshape  Development 


Graphs  of  Lineshapes 


Cross  Section  Development  - >  Graphs  of  Cross  Sections 


\ 


Overlay  Line  shape/ Pump  Graph 


Rate  Equation  Analysis 


ODE  Solution  Routine 


Solution  Selection  Algorithm 


Interpolation, I nte gration^ and  Tabular  Routines 


>1  Output  Plotting  Routines  and  Graphs 


Figure  8.  The  architecture  used  to  develop  the  simulation  of  the  DPAL  model.  Oc¬ 
tagons  are  output  plots  or  printed  output  system  characteristics.  Rectangles  are  algo¬ 
rithms. 
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quency,  as  this  is  the  frequency  of  the  laser’s  output.  7 „(z)  is  plotted  for  the  system. 
The  average  value  of  7^  over  z  is  compared  to  a  and  should  be  found  to  be  roughly 
equivalent.  A  determination  of  the  degree  to  which  the  chosen  solution  matches  the 
boundary  conditions  is  provided. 
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V.  Results  and  Simulation  Comparisons 


5.1  Chapter  Overview 

To  validate  the  model  this  thesis  develops,  which  is  henceforth  known  as  the  three 
level  numerical  model;  it  must  be  tested  against  similar  previously  vetted  systems 
like  the  Lewis  model  and  the  Hager  model.  If  the  outputs  of  the  model  developed  are 
comparable  to  those  developed  by  Lewis  and  Hager  in  the  regime  where  those  models 
are  known  to  operate  well,  then  the  model’s  output  outside  of  those  regimes  is  more 
believable.  The  model  will  also  be  tested  against  experimental  results.  Further, 
the  model  will  be  shown  to  simulate  a  DPAL  with  an  extremely  high  initial  pump 
intensity  and  other  features  typically  found  only  during  pulsed  operation.  The  model 
will  also  perform  Rigord  analysis  for  a  general  DPAL  case.  Together  these  outputs 
will  show  the  effectiveness  and  utility  of  the  model  and  its  simulation. 

5.2  Comparison  to  Lewis  Model 

5.2.1  Inputs. 

The  inputs  of  Table  6  are  derived  from  the  DPAL  quasi  two  level  regime  of  Lewis’s 
thesis  [5] .  The  inputs  arc  listed  in  the  units  used  by  Lewis  in  [5]  rather  than  units  used 
by  the  three  level  numerical  model.  If  an  input  was  not  listed  by  Lewis,  a  suitable 
value  was  devised.  This  applies  specifically  to  the  temperature  of  the  alkali  for  plots 
provided  in  [5].  The  alkali  used  for  this  comparison  was  Rb,  which  is  used  throughout 
the  comparison  between  the  three  level  numerical  model  and  the  Lewis  model.  These 
inputs  will  be  used  throughout  the  comparison  to  Lewis’s  model  and  the  three  level 
numerical  model  for  the  remainder  of  this  section  unless  otherwise  noted.  Lewis  lists 
one  of  his  buffer  gases  as  ethane  in  [5].  However,  his  theortical  development  and  plots 
list  methane  as  the  buffer  gas  used.  Further  [5]  lists  methane,  not  ethane,  as  his 
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buffer  gas  on  all  other  tables.  Therefore,  methane,  not  ethane  will  be  used  for  this 
development  [5]. 


Table  6.  Simulation  Input  Parameters  from  Lewis  Model 


Input  Parameter 

Units  required 

Value 

Symbol 

Temperature  of  Cell 

K 

455 

T 

Temperature  of  Methane 

K 

455 

T1 

meth 

Temperature  of  Helium 

K 

455 

Uffe 

Temperature  of  Alkali 

K 

455 

Talk 

Partial  Pressure  of  Methane 

T  orr 

600 

Mmeth 

Partial  Pressure  of  Helium 

T  orr 

200 

Mue 

Partial  Pressure  of  Alkali 

T  orr 

0 

Maik 

Total  Alkali  Number  Density 

cm~6 

6.1  x  1012 

Nt 

Length  of  gain  medium 

cm 

8.00 

k 

Transmission  of  windows  around  gain 

Unitless 

1.0 

Tg 

High  Reflector  Reflectivity 

Unitless 

1.0 

Ri 

Output  Coupler  Reflectivity 

Unitless 

1.0 

R‘2 

Initial  Pump  Intensity  at  linccenter 

Wcm~2 

200 

IpO 

Diode  Pump  Line  Center  frequency 

THz 

384.23 

Vd 

Diode  Pump  FWHM  frequency 

GHz 

3 

^ Pfwhm 

Space  Between  Cavity  Mirrors 

cm 

50 

^mirror 

5.2.2  Cross  Section  Broadening  Comparison. 

Lewis  details  that  the  main  mechanism  for  broadening  in  the  DPAL  system  is  the 
pressure  of  the  buffer  gas  [5] .  Lewis  then  investigates  the  effects  of  pressure  broadening 
at  several  different  pressures  on  the  emission  cross  section  of  the  3  to  1  transition, 
03i  shown  in  Figures  9  and  10.  The  analogous  results  for  the  three  level  numerical 
model  are  provided  by  accepting  that  the  Lewis  model  neglects  the  terms  fiso  and 
fji  by  assuming  only  one  isotope  is  present  and  the  hyperhne  states  all  are  equally 
populated  in  Figures  11  and  12.  The  Lewis  model  uses  units  of  cm 2  as  opposed  to  m 2 
which  will  induce  a  shift  of  10-4  when  comparing  cross  sections  between  the  Lewis 
model  to  the  three  level  numerical  model.  Notice  that  Figures  9  and  11  differ  slightly 
most  likely  clue  to  a  difference  in  temperature  as  Lewis  did  not  list  the  temperature 
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he  used  to  create  Figure  9.  Figures  10  and  12  are  identical.  This  implies  that  the 
lineshapes  and  cross  sections  developed  by  the  three  level  numerical  model  agree  well 
with  those  of  the  Lewis  model  especially  at  high  pressures.  Note  Figures  9,  10,  11, 
and  12  are  all  offset  in  frequency  space  by  za.  i  . 


o"3i  Cross  Section 

(T)i  (cm2) 


Figure  9.  a 31  for  100  Torr  helium  and  100  Torr  methane  from  the  Lewis  model  [5] 


5.2.3  Absorption  Profile  Comparison. 

The  absorption  profiles  for  the  DPAL  regime  from  the  Lewis  model  are  given  by 
Figure  13  from  [5].  The  same  inputs  were  provided  to  the  three  level  numerical  model 
and  the  comparable  output  is  given  by  Figure  14.  Notice  that  the  frequency  axis  in 
Figure  13  is  offset  based  on  pump  line  center  frequency  z/d,  the  intensity  is  in  Wcm~2, 
and  the  position  within  the  gain  medium  (z)  is  in  cm.  Figure  13’s  units  are  based 
on  the  MKS  system.  Figure  14  shows  similar  features  to  Figure  13,  however,  small 
changes  from  the  more  complete  rate  equation  analysis  and  the  use  of  the  terms  fiSO 
and  fji  are  noticeable  at  the  far  end  of  the  cell.  Hence,  the  Lewis  model  captures  a 


36 


(T 31  Cross  Section 

cr,,  (cm2) 


Figure  10.  (731  for  1000  Torr  helium  and  1000  Torr  methane  from  the  Lewis  model  [5] 
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Figure  11.  <731  for  100  Torr  helium  and  100  Torr  methane  from  the  three  level  numerical 
model 
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Figure  12.  0-31  for  1000  Torr  helium  and  1000  Torr  methane  from  the  three  level 

numerical  model 


great  deal  of  the  absorption  effects  within  the  gain  medium,  but  it  does  not  develop 
a  wholly  accurate  picture  of  the  attenuation  of  the  pump  wave. 

The  Lewis  model  also  develops  the  attenuation  of  the  pump  intensity  under  the 
quasi-two  level  approach  (QTLA),  and  based  upon  the  QTLA  assumption,  is  able 
to  develop  the  amount  of  attenuation  of  the  pump  intensity  due  to  lasing.  This 
development  is  somewhat  ad  —  hoc ,  and  will  not  work  properly  at  threshold  and  will 
not  give  a  spectral  profile  for  the  occurrence  of  lasing.  That  is,  the  pump  is  assumed 
to  cause  lasing  to  occur  and  then  will  be  attenuated  to  a  greater  degree  based  upon 
an  assumed  lasing  intensity  which  will  occur  over  all  pump  frequencies  to  an  equal 
degree.  This  effect  can  be  seen  in  Figure  15  with  the  inputs  of  Table  6.  In  actuality 
only  those  pump  photons  which  cause  lasing  to  occur  in  the  gain  medium  will  see 
this  effect.  Hence,  only  certain  pump  frequencies  will  exhibit  this  effect;  those  that 
induce  lasing  to  occur.  Those  pump  frequencies  which  do  not  provide  enough  energy 
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Figure  13.  The  Lewis  model  3D  absorption  profile  without  lasing  for  inputs  of  Table 
6.  Note  the  units  of  the  plot  are  not  MKS  and  the  frequency  is  offset  by  vd  [5]. 


Figure  14.  The  three  level  numerical  model  3D  absorption  profile  without  lasing  for 
inputs  of  Table  6 


to  maintain  gain  above  loss  will  not  observe  this  effect.  This  can  be  seen  in  Figure 
16,  created  with  the  inputs  of  Table  6  by  the  three  level  numerical  model.  Lasing 
is  observed  to  occur  in  the  frequency  domain  between  the  two  peak  features  which 
propagate  through  to  the  end  of  the  cell.  Also,  note  that  it  appears  that  the  pump 
intensity  in  the  area  without  lasing  has  grown  between  Figures  14  and  16.  This 
is  simply  an  optical  illusion  due  to  the  degree  to  which  a  discontinuity  appears  in 
frequency  space  due  to  achieving  threshold  inversion.  A  close  inspection  of  Figures 
14  and  16  will  reveal  this  fact.  By  comparing  Figures  15  and  16,  one  can  see  that  the 
Lewis  model  is  able  to  only  approximate  the  effect  of  the  attenuation  of  the  pump  due 
to  intra-cavity  lasing.  This  effect  is  negligible  well  above  threshold  (30  times  Isat)- 
The  three  level  numerical  model  also  simulated  7(2)  for  the  lasing  region  of  Figure 
16  given  in  Figure  17.  Figure  17  only  shows  the  area  with  positive  gain.  An  effective 
laser  should  end  when  gain  dips  below  zero  as  the  intra-cavity  lasing  waves  will  be 
absorbed  beyond  this  point. 

5.3  Comparison  to  Hager  Model 

5.3.1  Spectral  Profile  Comparison. 

In  [3],  Hager  develops  a  spectral  line  profile  of  the  absorption  of  the  D1  transition 
for  all  hyperfine  states  of  Rb.  The  data  Hager  presents  is  based  upon  experimental 
data  and  is  then  fit  to  his  development  of  the  Voigt  profile.  The  absorption  profile 
calculated  by  Hager  is  given  in  Figure  18.  The  analogous  emission  spectra  is  provided 
for  the  three  level  numerical  model  in  Figure  19.  Figure  18’s  absorption  features  are 
shown  to  be  exactly  mirrored  in  the  emission  profile  given  by  Figure  19.  Thus,  the 
three  level  model  is  able  to  fit  both  experimental  lineshape  data  and  the  Hager  model 
lineshape  for  a  spectra  including  all  of  the  naturally  occurring  isotopes  of  Rb  and  all 
of  the  hyperfine  transitions  of  the  D1  manifold. 


41 


3D  Absorption  F’rofilc  (QTLA) 

0  >  cm 

4 

>v  6 


150 

100  W 

I( - > 

50  cm* 

0 


Figure  15.  The  Lewis  model  3D  absorption  profile  with  QTLA  lasing  inputs  of  Table 
6.  Note  the  units  of  the  plot  are  not  MKS  and  the  frequency  is  offset  by  vd  [5]. 
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Figure  16.  The  three  level  numerical  model  3D  absorption  profile  with  lasing  for  inputs 
of  Table  6.  Lasing  is  only  occurring,  in  frequency  space,  in  the  region  between  the  two 
peak  features  which  propagate  throughout  the  cell. 
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Figure  17.  The  three  level  numerical  model  determination  of  7 (z)  for  inputs  of  Table  6 
within  the  lasing  region.  Notice  that  the  gain  is  only  provided  while  7  is  above  zero. 
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Figure  18.  The  hyperflne  absorption  profile  for  the  Rb.  D1  manifold  offset  by  z^i  [3] 
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Figure  19.  The  three  level  numerical  model  lineshape  for  the  Rb.  D1  manifold  offset 

by  v2i 
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5.4  CW  Simulation  of  a  Pulsed  System 


Pulsed  DPAL  systems  typically  operate  at  extremely  high  intensities  compared  to 
CW  systems,  however,  as  long  as  the  population  densities  reach  their  equilibrium  val¬ 
ues  a  CW  simulation  is  apt  for  a  pulsed  system.  Even  in  systems  which  do  not  achieve 
equilibrium,  the  rate  equations  for  the  population  concentrations  remain  unchanged 
from  Chapter  II.  So,  if  the  rate  at  which  populations  change  with  time  is  small  with 
respect  to  the  time  scale  of  the  pulse  width  of  the  diode  and  the  populations  are 
assumed  to  reach  semi-equilibrium  quickly  after  interaction  with  the  pulse,  then,  this 
development  is  still  at  least  somewhat  valid.  Hence,  the  three  level  numerical  model 
can  be  applied  to  high  intensity  systems  and  may  be  used  to  give  rough  estimates  of 
the  characteristics  of  some  pulsed  DPAL  systems.  Typical  inputs  for  the  operation 
of  a  pulsed  DPAL  system  can  be  found  in  Table  7. 


Table  7.  Simulation  Input  Parameters  for  CW  Simulation  of  a  Pulsed  System 


Input  Parameter 

LInits  required 

Value 

Symbol 

Temperature  of  Cell 

K 

500 

T 

Temperature  of  Methane 

K 

500 

nr 

- L  meth 

Temperature  of  Helium 

K 

500 

THe 

Temperature  of  Alkali 

K 

500 

Talk 

Partial  Pressure  of  Methane 

T  orr 

1000 

Mmeth 

Partial  Pressure  of  Helium 

T  orr 

1000 

HI  He 

Partial  Pressure  of  Alkali 

T  orr 

0 

Malk 

Total  Alkali  Number  Density 

m 

3.79  x  1019 

Nt 

Length  of  gain  medium 

m 

0.01 

h 

Transmission  of  windows  around  gain 

Llnitlcss 

1.0 

t9 

High  Reflector  Reflectivity 

Llnitlcss 

1.0 

R\ 

Output  Coupler  Reflectivity 

Unitlcss 

0.5 

Rz 

Initial  Pump  Intensity 

GWm~2 

5 

IpO 

Diode  Pump  Line  Center  frequency 

THz 

384.23 

Vd 

Diode  Pump  FWHM  frequency 

GHz 

50 

^ Pfwhm 

Space  Between  Cavity  Mirrors 

m 

0.1 

^mirror 

The  outputs  of  the  three  level  numerical  model  for  the  inputs  of  Table  7  are 
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given  in  Figures  20  -  23  and  Table  8.  Figure  20  provides  the  degree  to  which  the 
diode  pump  lineshape  was  matched  with  the  lineshape  of  the  alkali  atom.  In  this 
case,  the  areas  overlap  about  80  percent.  In  general,  the  pump  lineshape  should  be 
well  matched  to  the  transition  it  is  attempting  to  pump,  otherwise  energy  will  be 
lost.  Based  on  Figure  21  most  of  the  pump’s  input  intensity  is  still  present  at  the 
output  coupler,  implying  that  the  gain  cell  should  be  extended  if  possible.  Unlike  in 
Figure  16,  the  pump  intensity  in  Figure  21  is  not  observed  to  have  a  region  which 
is  being  attenuated  by  laser  operation  and  a  region  which  is  not.  This  is  due  to 
lasing  occurring  across  the  entire  pump  spectrum  provided  in  Figure  21.  Figure  22 
shows  the  gain  as  a  function  of  z.  The  gain  can  be  seen  to  decrease  as  the  pump 
attenuates.  However,  the  gain  does  not  decrease  linearly  with  position,  which  implies 
that  approximations  such  as  LAND  may  not  be  valid  in  this  case.  The  total  gain 
in  Table  8  is  the  average  value  of  the  gain  in  Figure  22.  In  Table  8  7 (total) /a  is  a 
measure  to  how  well  the  system  followed  the  approximation  that  gain  equals  loss.  At 
infinite  fidelity  7 {total) /a  should  equal  unity.  So,  7 [total) /a  is  a  measure  of  merit 
of  how  well  the  system  performed.  The  average  residual  in  solutions  listed  in  Table 
8  is  the  of  average  of  all  residuals  of  the  routine  which  chose  the  best  adherence  to 
the  boundary  condition  at  the  output  coupler,  i.e.  Equation  (24).  This  average  is 
only  computed  for  those  cases  in  which  lasing  is  determined  by  the  simulation  to  be 
occurring.  Figure  23  shows  the  forward  traveling  wave  (upper)  and  the  backward 
traveling  wave  (lower).  Notice  that  difference  between  the  upper  wave  and  the  lower 
wave  at  the  end  of  the  cell  is  equal  to  the  reflectivity  of  the  output  coupler  times  the 
square  of  the  transmissivity  of  the  gain  cell,  which  in  this  case  is  0.5. 
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Figure  20.  The  lineshape  g$i  with  the  pump  lineshape  gp  overlayed  to  show  the  degree 
of  area  matching.  Note  the  frequency  is  offset  by  [5]. 


Pump  Intensity  vs.  Freq.  vs.  Pos. 
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Figure  21.  The  intensity  of  the  pump  (Ip)  as  it  propagates  through  the  cell 
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y  (Integrated)  vs.  Pos 


Figure  22.  The  gain  7  as  a  function  of  z. 


Figure  23.  The  Integrated  Plus  Wave  Intensity(upper)  and  Minus  Wave  Inten- 
sity(lower)  as  functions  of  z.  The  output  coupling  is  the  space  between  the  upper 
wave  and  lower  wave  at  the  edge  of  the  graph. 
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Table  8.  Simulation  Outputs  Characteristics  for  CW  Simulation  of  a  Pulsed  System 


Input  Parameter 

Units 

Value 

Symbol 

Output  Laser  Intensity 

MWrrr 2 

63.7 

lout 

Output  Laser  Frequency 

THz 

377.107 

n 

lout/  Ip 

Out  of  1 

0.128 

n/a 

Ip/ 1 sat 

n/a 

5002. 

n/a 

7  (total) /a 

n/a 

1.001 

n/a 

Average  residual  in  solutions 

n/a 

0.0094 

n/a 

5.5  Simulation  of  a  System  Near  Threshold 

One  of  the  more  difficult  regimes  to  model  in  most  cases  is  the  regime  at  or  near 
threshold.  The  three  level  numerical  model  is  able  to  simulate  threshold  systems 
effectively.  The  inputs  for  three  level  numerical  model  for  a  threshold  system  are 
given  by  Table  9.  Under  these  conditions,  only  a  the  center  bandwidth  of  the  pump 
wave  can  induce  a  population  inversion,  and  thus  lasing.  This  effect  is  visible  in  Figure 
24.  The  outputs  of  the  three  level  numerical  model  are  provided  in  Table  10.  Near 
threshold  the  three  level  numerical  model  had  difficulty  obtaining  an  interpolation 
for  7 {total)  so  it  is  not  listed.  The  difficulty  seems  to  arise  from  the  sharpness  of  the 
gain  profile  for  a  threshold  case. 
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Table  9.  Simulation  Input  Parameters  for  a  Threshold  System 


Input  Parameter 

Units  required 

Value 

Symbol 

Temperature  of  Cell 

K 

450 

T 

Temperature  of  Methane 

K 

450 

Tmeth 

Temperature  of  Helium 

K 

450 

THe 

Temperature  of  Alkali 

K 

450 

Talk 

Partial  Pressure  of  Methane 

T  orr 

500 

Mmeth 

Partial  Pressure  of  Helium 

T  orr 

500 

MHe 

Partial  Pressure  of  Alkali 

T  orr 

0 

Malk 

Total  Alkali  Number  Density 

m 

6.1  x  1018 

Nt 

Length  of  gain  medium 

m 

0.02 

l9 

Transmission  of  windows  around  gain 

Unitless 

1.0 

Tg 

High  Reflector  Reflectivity 

Unitless 

1.0 

R\ 

Output  Coupler  Reflectivity 

Unitless 

0.9999 

R‘2 

Initial  Pump  Intensity 

kWm~'2 

600 

IP0 

Diode  Pump  Line  Center  frequency 

THz 

384.23 

Vd 

Diode  Pump  FWHM  frequency 

GHz 

5 

^ Pfwhm 

Space  Between  Cavity  Mirrors 

m 

0.1 

dmirror 

Table  10.  Simulation  Outputs  Characteristics  for  a  Threshold  System 


Input  Parameter 

LInits 

Value 

Symbol 

Output  Laser  Intensity 

Wm~ 2 

11.2 

d-out 

Output  Laser  Frequency 

THz 

377.109 

O 

d out  /  dp 

Out  of  1 

0.000020 

n/a 

dp/ 1 sat 

n/a 

1.25 

n/a 

'y  (total)/ a 

n/a 

n/a 

n/a 

Average  residual  in  solutions 

n/a 

0.35 

n/a 
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Figure  24.  The  laser  output  intensity’s  spectral  width  within  the  pump  profile.  To  get 
total  output  intensity  one  must  multiply  by  gp  and  integrate  over  all  space. 
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VI.  Conclusions 


6.1  Comparison  to  Other  Models 

The  three  level  numerical  model  is  able  to  replicate  the  results  of  both  the  Lewis 
model  and  the  Hager  model  and  to  produce  expected  results  that  neither  of  these 
models  are  able  to.  The  three  level  numerical  model  offers  a  much  higher  fidelity  under 
many  circumstances,  but  with  that  added  fidelity  comes  a  much  more  cumbersome 
development  and  a  lack  of  intuitive  understanding  of  the  problem,  which  is  a  hallmark 
of  the  Lewis  and  Hager  models.  Hence,  the  three  level  numerical  model  provides 
another  option  for  the  simulation  of  DPAL  systems,  but  does  not  supersede  previous 
developments. 

6.2  Use  as  a  Research  Tool 

Though  the  three  level  numerical  model  can  be  somewhat  difficult  to  utilize  and 
to  interpret,  it  provides  a  great  deal  of  fidelity  for  research  into  many  areas  which  no 
other  DPAL  simulation  can  provide.  The  three  level  numerical  model  handles  broad¬ 
band  pumping,  produces  a  complete  analysis  of  the  rate  equations  for  the  three  level 
DPAL  system,  selects  preferential  gain  from  cavity  mode  spacing,  develops  the  hy- 
perfine  lineshape  for  multiple  isotopes,  allows  the  use  of  multiple  buffer  gases,  allows 
for  quenching  to  be  simulated,  and  solves  a  set  of  three  coupled  non-linear  transcen¬ 
dental  differential  equations  to  characterize  any  three  level  CW  DPAL  system  under 
any  regime  to  high  fidelity.  Hence,  the  three  level  numerical  model  is  an  excellent 
tool  for  the  simulation  and  characterization  of  any  CW  DPAL  system.  Further,  the 
three  level  numerical  model  can  be  used  as  a  tool  to  approximate  the  parameters  of 
many  pulsed  systems. 
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6.3  Future  Model  Development 


6.3.1  Mode  Volume. 

The  three  level  numerical  model  does  not  account  for  the  incomplete  filling  of 
cavity  mode  volumes  by  the  pump  and  currently  does  not  have  the  ability  to  simulate 
the  lost  energy  from  this  effect.  Any  further  iteration  of  the  three  level  numerical 
model  should  include  this  effect  as  it  is  a  persistent  issue  for  DPAL  systems  and 
without  its  effects  a  model  cannot  hope  to  completely  capture  the  effects  observed 
experimentally.  The  main  hurdle  to  implementation  is  the  lack  of  references  on  the 
subject  of  mode  volume  characterization  for  DPAL  systems. 

6.3.2  Pulsed  DPAL  Systems. 

Probably  the  most  important  addition  to  the  three  level  numerical  model  would 
be  the  addition  of  time  dependence.  Though  this  would  add  a  great  deal  of  com¬ 
plexity  any  high  power  system  will  most  likely  utilize  pulsed  operation.  Much  of  the 
current  research  in  DPAL  systems  involves  the  use  of  pulse  operated  DPAL  system. 
Therefore,  any  further  development  of  the  model  should  include  the  ability  for  pulsed 
operation.  While  theoretical  models  do  exist  for  the  pulsed  operation  of  lasers  the 
inherent  complexity  of  a  time  dependent  laser  system  forces  super-computing  as  a 
near  necessity  for  any  high  fidelity  systems. 
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Appendix  A. 


Mathematica  Code  to  Solve  Rate  Equations 
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DPAL  3-Level  Modeling  Rate  Equation  Solution 


I.  Rate  Equations 


ClearAll[Nl,  N2,  N3,  Reqnl,  Reqn2 ,  Reqn3,  Reqn4 ,  f,  z] 

Reqnl  [Nl_,  N2_,  N3_]  :  =  N3  *  A31  [iso]  +  N2  *  A21  [iso]  +  N3  *  k31  +  N2  *  k21  + 

Ip[z]  * glpOloop  Ip[z]  * glpOloop 

-B13  [vpump  ,  iso]  *  N1  *  -  +  B31  [vpump  ,  iso]  *  N3  *  — 


B21[vlaser,  iso]  * 


c  /  ng 

Il[z]  * gapproxfixtota!211oop 


c  /  ng 


c  /  ng 


g2 

N2 - N1 

gi 


Reqn2  [Nl_,  N2_,  N3_]  :=  N3  *  k32  [ species ,  iso]  *M[species]  -  N2  *  k23  [ species ,  iso]  *M[species] 


N2  *  A21  [iso]  -  N2  *  k21  -  B21  [vlaser  ,  iso]  * 


II  [z]  * gapproxfixtota!211oop 


c  /  ng 


g2 

N2 - N1 

gi 


Reqn3  [Nl_,  N2_,  N3_]  :  = 

-N3 * k32 [species ,  iso]  *M[species]  +  N2  *  k23  [ species ,  iso]  *M[species]  -N3*A31[iso]  - 

Ip[z]  * glpOloop  Ip[z]  * glpOloop 

N3  *  k31  +  B13  [vpump  ,  iso]  *  -  *  N1  -  B31  [vpump  ,  iso]  *  N3  *  - 

c  /  ng  c  /  ng 


Reqn4  [Nl_,  N2_,  N3_]  :  =  N1  +  N2  +  N3 
N2 

f  =  - ; 

N2  +  N3 

Solve  [  {Reqnl  [Nl,  N2 ,  N3]  ,  Reqn2[Nl,  N2 ,  N3]  ,  Reqn4[Nl,  N2,  N3]  }  ==  {0,  0,  Nt},  {Nl,  N2,  N3}] 


A/1  -»  A/t  +  (g2  gapproxf  ixtotal211oop  ngA/tS21  [vlaser,  iso]  J1  [  z  ]  )  / 

/  gapproxfixtotal211oop  ng  S21  [vlaser ,  iso]  Jl[z] 

cgl  -  k2 1  -5121  [iso]  -  - 


g2  gapproxf ixtotal211oop  ng$21  [vlaser,  iso]  J1  [z] 


c  gl 


-  k23 [ species  ,  iso]  A1[species] 
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Printed  by  Mathematicafor  Students 
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c  gl 


g2  gapproxf ixtotal2 Hoop  ng  A/t  £21  [vlaser,  iso]  J1  [z] 


k.2 1  +  Ji2 1  [  iso  ]  + 


gapproxf ixtota!2 lloop  ng£21  [vlaser,  iso]  II  [z] 


g2  gapproxf ixtotal211oop  ng£21  [vlaser,  iso]  II  [z] 

c  gl 

gJpOloop  ng  £13  [  vpump ,  iso]  Jp[z]\ 


A/t  -  - 


g2  gapproxf ixtotal211oop  ng£21  [vlaser,  iso]  J1  [z] 

c  gl 


gJpOloop  ng£13  [vpump,  iso]  Jp  [z] 


-  k2 1  -^?21  [iso]  - 


gapproxf ixtotal211oop  ng£21  [vlaser,  iso]  J1  [z] 


g2  gapproxf ixtotal2 lloop  ng£21  [vlaser,  iso]  J1  [z] 

c  gl 

k23[species,  iso]  M [ species ] j j j 

(  g2  gapproxf ixtotal2 lloop  ng£21  [vlaser,  iso]  J1  [z] 

k31  +  H31  [iso]  +  -  + 

\  c  gl 

gJpOloop  ng  £13  [vpump  ,  iso]  Jp  [  z  ]  gJpOloop  ng  £31  [  vpump  ,  iso]  Jp[z] 


-k21  -121  [iso]  - 


c  c 

gapproxf ixtotal2 lloop  ng£21  [vlaser,  iso]  J1  [z] 


g2  gapproxf ixtotal2 lloop  ng£21  [vlaser,  iso]  J1  [z] 

c  gl 


•  k23 [ species ,  iso]  A1[species] 


gapproxf ixtotal211oop  ng£21  [vlaser,  iso]  J1  [z] 

k21  +  K21  [iso]  +  - 

c 

g2  gapproxf ixtotal211oop  ng£21  [vlaser,  iso]  J1  [z] 

-  + 

c  gl 

gJpOloop  ng£13  [vpump,  iso]  Jp  [z]  \ 


g2  gapproxf ixtotal211oop  ng£21  [vlaser,  iso]  J1  [z] 

c  gl 
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Printed  by  Mathematicafor  Students 
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k32 [species,  iso]  M [ species ] J J  - 
[  g2  gapproxf ixtota!211oop  ngMl  [vlaser,  iso]  II  [z] 


c  gl 


+ k32 [ species ,  iso]  A1[species] 


1 


c  gl 


g2  gapproxf  ixtotal211oop  ng  Ml  Ml  [vlaser,  iso]  II  [z] 


gapproxf ixtotal211oop  ngMl  [vlaser ,  iso]  II  [z] 


k21  +  Ml  [iso]  + 


g2  gapproxf ixtotal211oop  ngMl  [vlaser,  iso]  II  [z] 


Nt 


c  gl 

gJpOloop  ng  M3  [  vpump  ,  iso]  Jp[z] 


/  g2  gapproxf ixtotal211oop  ngMl  [vlaser,  iso]  II  [z] 

V  c  gl 

gJpOloop  ng  M3  [  vpump  ,  iso]  Jp[z]\ 


(  gapproxf ixtotal211oop  ngMl  [vlaser,  iso]  II  [z] 

-k21-^21[iso]  -  - 

V  c 

g2  gapproxf ixtotal211oop  ngMl  [vlaser,  iso]  II  [z] 

c  gl 

k23 [species,  iso]  yM[species] 

[  gapproxf ixtotal211oop  ngMl  [vlaser,  iso]  II  [z] 

- k21  -  Ml  [iso] - 


g2  gapproxf ixtotal2 Hoop  ngMl  [vlaser,  iso]  II  [z] 

c  gl 


-  k23 [ species ,  iso]  A1[species] 


k31  +  30)1  [  iso  ]  + 


g2  gapproxf ixtotal211oop  ngMl  [vlaser,  iso]  II  [z] 

c  gl 


gJpOloop  ng  M3  [  vpump  ,  iso]  Jp[z]  gJpOloop  ng  Ml  [ vpump  ,  iso]  Jp[z] 

-  +  - 

c  c 

gapproxf ixtotal211oop  ngMl  [vlaser,  iso]  II  [z] 

-k21  -  Ml  [iso]  - - 


g2  gapproxf ixtotal211oop  ngMl  [vlaser,  iso]  II  [z] 


c  gl 


-  k23 [ species ,  iso] 


58 


Printed  by  Mathematicafor  Students 


4  |  Rate  Eqns  Fapp.nb 


M [ species ] 


k2 1  +  X2 1  [  iso  ]  + 


gapproxfixtota!211oop  ng  $21  [ vlaser  ,  iso]  XI  [z] 


g2  gapproxf ixtotal211oop  ng$21 [vlaser,  iso]  XI  [z] 

c  gl 

g/pOloop  ng  SI 3  [ vpump  ,  iso]Xp[z]\ 


g2  gapproxfixtotal211oop  ng$21  [vlaser,  iso]  XI  [z] 

c  gl 

k32 [species,  iso]  A1[species] 

Afi  ->  -  (g2  gapproxf  ixtotal 2 Hoop  ng  A/t  $21  [  vlaser  ,  iso]  XI  [  z]  )  j 

I  gapproxf ixtotal211oop  ng$21 [vlaser,  iso]  XI [z] 

cgl  -k21-.X21[iso]  -  - 


g2  gapproxfixtotal211oop  ng  $21  [vlaser  ,  iso]  XI  [z] 

c  gl 

k23 [ species ,  iso]  M [ species ] j  j  + 
g2  gapproxf ixtotal211oop  ngS21  [vlaser,  iso]  XI  [z] 


c  gl 


+  k32 [ species  ,  iso]  A1[species ] 


1 


c  gl 


g2  gapproxf ixtotal211oop  ng  Nt  S21  [vlaser,  iso]  XI  [z] 


gapproxf ixtotal211oop  ng  S21 [vlaser ,  iso]  XI  [z] 


k21+X21[iso]  + 

c 

g2  gapproxf ixtotal211oop  ng$21  [vlaser,  iso]  XI  [z] 

c  gl 

gXpOloop  ng  $13  [vpump ,  iso]Xp[z]\ 


At 


j  g2  gapproxf ixtotal211oop  ng$21  [vlaser,  iso]  XI  [z] 

c  gl 

gXpOloop  ng  $13  [vpump ,  iso]Xp[z]\ 


-k21-$21[iso]  - 


gapproxfixtotal211oop  ng  $21  [vlaser ,  iso]  XI [z] 
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g2  gapproxf ixtotal211oop  ngMl  [vlaser ,  iso]  II  [z] 

c  gl 


k23  [species,  iso]  yvt[species] 


gapproxf ixtotal211oop  ngMl  [vlaser,  iso]  II  [z] 

-k21  -  ^2 1  [iso]  - - 

c 

g2  gapproxf ixtotal2 Hoop  ngMl  [vlaser,  iso]  II  [z] 

- k23 [species,  iso]  A1[species] 

c  gl 

g2  gapproxf ixtotal211oop  ngMl  [vlaser,  iso]  II  [z] 

k31  +  H31  [iso]  +  - + 

c  gl 

gJpOloop  ng  M3  [ vpump  ,  iso]  Jp[z]  gJpOloop  ng  Ml  [vpump  ,  iso]  Jp[z] 

-  +  - 

c  c 

gapproxf ixtotal2 lloop  ngS21  [vlaser,  iso]  II  [z] 

-k21-Ml[iso]  - - 

c 

g2  gapproxf ixtotal211oop  ngMl  [vlaser,  iso]  II  [z] 

- k23 [species,  iso] 

c  gl 

'  gapproxf ixtotal211oop  ngMl  [vlaser,  iso]  II  [z] 

A1[species]  -  k21  +  M 1  [iso]  +  -  + 

,  c 

g2  gapproxf ixtotal211oop  ngS21  [vlaser,  iso]  II  [z] 

c  gl 

gJpOloop  ng  M3  [ vpump  ,  iso]  Jp[z] 
c 

g2  gapproxf ixtotal211oop  ngS21  [vlaser,  iso]  II  [z] 

c  gl 


k32 [species,  iso]  Mfspecies] 

I  1 

N3  ->  - - g2  gapproxf ixtotal2 lloop  ng  /Vt  Ml  [vlaser  ,  iso]  II  [z] 

V  c  gl 

gapproxf ixtotal2 lloop  ngMl  [vlaser,  iso]  II  [z] 

k21  +  K2 1  [iso]  +  -  + 

c 

g2  gapproxf ixtotal211oop  ngMl  [vlaser,  iso]  II  [z] 

c  gl 

gJpOloop  ng  M3  [  vpump  ,  iso]  Jp[z] 
c 
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Nt 


j  g2  gapproxf ixtotal211oop  ngi821  [vlaser,  iso]  II  [z] 

c  gl 

gJpOloop  ng  SI 3  [ vpump  ,  iso]  Jp[z]\ 


(  gapproxf ixtotal2 Hoop  ngi321  [vlaser,  iso]  II  [z] 

-k21  -  H21  [iso]  - - 

[  c 

g2  gapproxf ixtotal2 lloop  ngS21  [vlaser,  iso]  II  [z] 

c  gl 

k23 [species,  iso]  M [ species ] j j j 

f  g2  gapproxf ixtotal211oop  ng£21  [vlaser,  iso]  II  [z] 

k31  +  H31  [iso]  +  -  + 

\  c  gl 

gJpOloop  ng  S13  [ vpump  ,  iso]  Jp[z]  gJpOloop  ng  S3 1  [vpump  ,  iso]  Jp[z] 


c  c 

gapproxf ixtotal2 lloop  ngi321  [vlaser,  iso]  II  [z] 


-  k2 1  -H21  [iso]  - 


g2  gapproxf ixtotal2 lloop  ng.S21  [vlaser,  iso]  II  [z] 


c  gl 

gapproxf ixtota!211oop  ng£21  [vlaser,  iso]  II  [z] 


•  k23 [ species ,  iso]  A1[species] 


k21  +  K2 1  [iso]  + 


g2  gapproxf ixtotal2 lloop  ng^21  [vlaser,  iso]  II  [z] 


c  gl 

gJpOloop  ngS13  [vpump,  iso]  Jp  [z] 


g2  gapproxf ixtotal2 lloop  ngS21  [vlaser,  iso]  J1  [z] 

c  gl 

k32 [species,  iso]  M [ species ] j j jj 


II.  Full  Simplify  of  f 

f  is  found  to  be  a  function  of  M,  v,  and  z 
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N3  :  = 


-  B13  Ip 


B21  ^  II 

-  k21  -  A21  -  B21  II - k23Eth  MEth 


Nt  + 


B21  £  II  (-B13  Ip  +  k23Eth  MEth)  Nt 


/ 


B21 ^  II 

-  k21  -  A21  -  B21  II - k23Eth  MEth 


( -  k31  -  A31  -  B13  Ip  -  B31  Ip  -  k32Eth  MEth)  + 


B21 ^  II 


+  k32Eth  MEth 


N2  :  : 


( -  B13  Ip  +  k23Eth  MEth) 

B13  Ip  Nt 

-B13  Ip  +  k23Eth  MEth 

B21 ^  II 

-B13  Ip  I  -k21  -  A21  -  B21  II - k23Eth  MEth 


+  ( -  k31  -  A31  -  B13  Ip  -  B31  Ip  -  k32Eth  MEth) 


B21  ^  II  (-B13  Ip  +  k23Eth  MEth)  Nt 


Nt  + 


)/ 


B21  ^  II 

-  k21  -  A21  -  B21  II - k23Eth  MEth 


(-B13  Ip  +  k23Eth  MEth) 

(-  k31  -  A31  -  B13  Ip  -  B31  Ip  -  k32Eth  MEth)  + 


(-B13  Ip  +  k23Eth  MEth) 


B21 ^  II 


+  k32Eth  MEth 


N1  :  =  -  ( -  k21  k31  Nt  -  k31  A21  #  Nt  -  k21  A31  Nt  -  A21  A31  Nt  -  k31  B21  II  Nt  - 

A31  B21  II  Nt  -  k21  B31  Ip  Nt  -  A21  B31  Ip  Nt  -  B21  B31  II  Ip  Nt  - 
k23Eth  k31  #  MEth  Nt  -  k21  k32Eth  MEth  Nt  -  k32Eth  A21  MEth  Nt  - 
k23Eth  A31  MEth  Nt  -  k32Eth  B21  #  II  MEth  Nt  -  k23Eth  B31  Ip  MEth  Nt)  / 

(k21  k31  +  k31  A21  #  +  k21  A31  +  A21  A31  +  k31  B21  II  +  A31  B21  II  + 

k31  B21  ^  I1+A31B21^  II  +  k21  B13  Ip  +  A21  B13  Ip  +  k21  B31  #  Ip  +  A21  B31  Ip + 

B13  B21  #  II  Ip  +  B21  B31  II  Ip  +  B21  B31  ^  II  Ip  +  k23Eth  k31  #  MEth  +  k21  k32Eth  MEth  + 
k32Eth  A21  MEth  +  k23Eth  A31  MEth  +  k32Eth  B21  II  MEth  +  k23Eth  B21  ^  II  MEth  + 
k32EthB21^  II  MEth  +  k23Eth  B13  Ip  MEth  +  k32Eth  B13  Ip  MEth  +  k23Eth  B31  Ip  MEth) 
FullSimplify [ 


f] 


(S21  g2  II  (k31  +5131  +  S31  Ip)  +  k32Eth  (S21  g2  II  +  S13  gl  Ip)  AEth  )  / 
(  (k21  +  521  )  SI  3  gl  Ip  +  S21  II  (S13  gl  Ip  +  g2  (k31  +  5131  +  S31  Ip))  + 
(k23Eth  +  k32Eth )  (S21  g2  II  +  S13  gl  Ip  )  AEth  ) 


III.  Limiting  Cases  for  f 


1 .  Case  Ip  approaches  0 

Limit  [f.  Ip  ->  0] 

k31  +  5131  +  k32Eth  AEth 
k31  +  5131  +  (k23Eth  +  k32Eth)  AEth 
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2.  Case  II  approaches  0 

Limit  [  f ,  1 1  -»  0  ] 

k32Eth  AEth 

k21  +521  +  (k23Eth  +  k32Eth)  AEth 


3.  Case  Ip  approaches  0  and  II  approaches  0  in  said  order 

FullSimplify  [Limit  [Limit  [f.  Ip  ->  0],  II  ->  0]] 

k31  +  5131  +  k32Eth  AEth 
k31  +  5131  +  (k23Eth  +  k32Eth)  AEth 

4.  Case  Ip  approaches  0  and  II  approaches  0  in  said  order 

FullSimplify  [Limit  [Limit  [f,  II  -»  0],  Ip  ->  0]] 

k32Eth  AEth 

k21  +521  +  (k23Eth  +  k32Eth)  AEth 


5.  Case  MEth  approaches  oo 

FullSimplify  [Limit  [f ,  MEth  -»  oo  ]  ] 

k32Eth 

k23Eth  +  k32Eth 


It  is  of  note  that  if  A31,  k31,  A21,  k21  are  relatively  small  (which  they  should  be)  and  Ip  and  II  are  zero  then  the  same  result  is 
obtained  for  all  5  cases. 
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Three  Level  DPAL  CW  Model  for  Rb 


I.  User  Inputs  (in  mks  units) 


Temperatureofcell  =  450 (*K  of  cell*) ; 

TemperatureMeth  =  450 (*K  of  relaxtion  gas  Methane*); 

TemperatureHe  =  450 (*K  of  relaxtion  gas  Helium*); 

Alkalitemperature  =  450 (*K*); 

TotalAlkaliConcentration  =  1.5  *  1019  (*irT3*)  ; 

Celllength  =0.03 (*m*) ; 

Celltransmission  =  1 . 0 (*Unitless*) ; 

HRref lectivty  =  1 . 0 (*Unitless*) ; 

OCref lectivty  =  0 . 99 (*Unitless*) ; 

IntialPumpIntesnity  =  50  000  000  (*In  w/m2*)  ; 

PartialPressureMethane  =  10  *  6666 (*Pa*) ; 

PartialPressureAlkali  =  0(*Pa*); 

PartialPressureHelium  =  10  *  6666  (*Pa*) ; 

Pumplinecenter  =  384.2304844685  *  1012  (*frequency  of  pump  line  center  in  Hz*); 

PumpFWHM  =  25  *  109(*the  FWHM  of  the  pump  in  Hz*) ; 

DistanceBetweenMirrors  =  0.1(*in  m*) ; 

Fidelityl  =  10 

(*The  Amount  of  Grid  Points  in  I  for  3-D  Grid  on  which  to  place  solution*) ; 

Fidelityv  =  10(*The  Amount  of  Grid  Points  in  v  pump  for  3- 
D  Grid  on  which  to  place  solution*) ; 

Fidelityz  =  10(*The  Amount  of  Grid  Points  in  z  for  3-D  Grid  on  which  to  place  solution*) 

II.  Constants 


A.  Common  Physical  Constants  (in  standard  SI  units  mks) 

c  =  299  792  458  (*Wikipedia  m/s*) ; 
kb  =  1.380650424  *  10~23  (*  J/K*)  ; 
h  =  6.6260689633  *  10'34  (* J/Hz*)  ; 
amu  =  1 . 6  6  0  5  3  8  7  8  2  8  3  *  10  27  (*kg*)  ; 

B.  Einstein  A  Coeffcients 


A21 [85] 

A21 [87] 

A31 [85] 

A31 [87] 

A32 [85] 
A32 [87] 


1 

-  (*Lewis  Hz*) ; 

27 . 679  *  10‘9 

1 

- (*Lewis  Hz*); 

27.7*  10~9 

1 

- (*Lewis  Hz*) ; 

26.2348  *  10“9 

1 

- (*Lewis  Hz*) ; 

26.24  *  10'9 

0; 

0; 
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III.  Parameters 

A.  Enviromental  Paremeters 

temp  =  294 . 15  (*K*) ; 

tempcell  =  Temperatureof cell ; 

celloutsidepressure  =  101.3  * 103  (*Pa*) ; 

B.  Design  Parameters 

celllength  =  Celllength (*m*) ; 

celltransmission  =  Celltransmission (*Unitless*) ; 
hrref lectivty  =  HRreflectivty (*Unitless*) ; 
ocref lectivty  =  OCreflectivty (*Unitless*) ; 

Atomicmass [ 87 ]  =  86.909180527  (*Lewis  amu*) ; 

Atomicmass [85]  =  84 . 911789738 (*Lewis  amu*); 

n=  1.  (*index  of  refraction  Unitless*) ; 

ng  =  1. (*  group  index  of  refraction  Unitless*) ; 

v21  =  377.1074635  * 1012  (*frequency  21  tranisition  Lewis  in  Hz*)  ; 
v31  =  384 .2304844685  *  1012  (*Lewis  Hz*); 
v32  =  7.123020968  * 1012(*Lewis  Hz*); 
distance  =  DistanceBetweenMirrors ; 

C.  Laser  and  Cell  Input  Parameters 

(*If  you  wish  to  add  more  species  you  must  do  so  in  II.  C.  ,  I.  D.,  and  in  V.  A.*) 
TotalAlkaliConcentration  =  TotalAlkaliConcentration  (*In  irT3*)  ; 

tempMeth  =  TemperatureMeth (*K*) ; 

MMeth  =  PartialPressureMethane  (*Partial  Pressure  of  Methane  in  Pa*)  ; 

MMeth  =  16.04246  *amu(*Mass  of  Methane  Lewis  kg*) ; 

■yMethl  = 218 304 . 576 (*For  Rb87  and  D1  manifold [ iso? ]  in  Hz/Pa  Hager*); 

■yMeth2  =  196 549 . 137 (*For  Rb87  and  D2  manifold [ iso? ] in  Hz/Pa  Hager*); 

<5Mel  =  -59459.4595;  (*Hz/Pa  Collision  induced  shift  D1  Hager*) 

<5Me2  =  -52552.5526;  (*Hz/Pa  Collision  induced  shift  D2  Hager*) 

tempAlk  =  Alkalitemperature  (*K*) ; 

MAlk  =  PartialPressureAlkali ;  (*Partial  Pressure  of  Alkali  in  Pa*) 

MAlk  =  85.4678;  (*Mass  of  Alkali  function  of  iso?  in  kg*) 

YAlkl  =  0(*  in  Hz/Pa*)  ; 

YAlk2  =  0(*  in  Hz/Pa*)  ; 

tempHe  =  TemperatureHe (*K*) ; 

MHe  =  PartialPressureHelium  (*Partial  Pressure  of  Helium  in  Pa*) ; 

MHe  =  4.002602  *amu(*Mass  of  Helium  in  kg*) ; 

■yHel  =  141 785 . 446 (*For  Rb87  and  D1  manifold [ iso? ]  in  Hz/Pa  Hager*); 

YHe2  =  150 037 . 509 (*For  Rb87  and  D2  manifold [ iso? ]  in  Hz/Pa  Hager*); 

<5Hel  =  35333.8335;  (*Hz/Pa  Collision  induced  shift  D1  Hager*) 

<5He2  =  2775.69392;  (*Hz/Pa  Collision  induced  shift  D2  Hager*) 
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D.  Diode  Input  Lineshape 


glpO  [v_]  :: 


7T  r  Pumplinecenter  ,  r  (v  -  Pumplinecenter ) 2  1 

—  PumpFWHM  1  +  Erf  -  Exp - ; 

2  V2~  PumpFWHM  2  PumpFWHM2 


(*This  is  the  Input  Lineshape  function  for 
gaussian  distribution  normalized  to  1  and  is  unitless*) 

(*glp0 [ v _ ] : =UnitStep [v-v31 +PumpFWHM/2] *UnitStep [v31 + PumpFWHM/ 2 -v] 

(*This  is  a  square  wave  and  may  be  substituted*)*) 

IpO  =  IntialPumpIntesnity  (*w/m2*)  ; 

Plot  [glpO  [v]  ,  {v,  v31  -  PumpFWHM  *  3,  v31  +  PumpFWHM  *  3}  , 

PlotLabel  -»  "gpump(v)  vs.  Freq .  ",  AxesLabel  ->  {gp[vp],  Hz}] 


IV.  Kinetics 


A.  k  Coeffcients  (quenching  rates  only) 

k21rate[l]  =0  /  (kb  * tempMeth)  (*  Methane  Hz  input  the  k  rate  in  m3/s*); 
k31rate[l]  =  0  /  (kb  *  tempMeth)  (*  Methane  Hz  input  the  k  rate  in  m3/s*)  ; 
k21rate[2]  =  0  /  (kb*tempHe)  (*  Helium  Hz  input  the  k  rate  in  m3/s*) ; 
k31rate[2]  =  0  /  (kb*tempHe)  (*  Helium  Hz  input  the  k  rate  in  m3/s*) ; 
k21rate[3]  =  0  /  (kb*tempAlk)  (*  Alkali  Hz  input  the  k  rate  in  m3/s*) ; 
k31rate[3]  =  0  /  (kb*tempAlk)  (*  Alkali  Hz  input  the  k  rate  in  m3/s*)  ; 
k21  :=  Sum [k2 Irate [ species ]  *M[species],  {species,  1,  speciesmax,  1}]; 
k31  :=  Sum [k31rate [ species ]  *M[species],  {species,  1,  speciesmax,  1}]; 

B.  Energy  Differences  D  manifold  states 

AE21  =  h  *  v21  (*J*)  ; 

AE31  =  h  *  v31  (*J*)  ; 

AE32  =  h  *  v32  (*J*)  ; 

C.  Boltzmann  Ratios 
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D.  k  Coeffcients  (for  partner  species  on  the  3  and  2  levels) 

(*If  you  wish  to  add  more  species  you  must  do  so  in  II.  C.  ,  I.  D.,  and  ine  V.  A.*) 
k32Meth[85]  =  3.16 *10'16/  (kb  * tempMeth)  (*Hager  provided  the  number  listed  in  mA3/s, 
but  the  value  of  the  coeffcient  is  in  Hz /Pa*) ; 
k32Meth[87]  =  3.16  *10"16/  (kb  * tempMeth)  (*Hager  provided  the  number  listed  in  mA3/s, 
but  the  value  of  the  coeffcient  is  in  Hz /Pa*) ; 
k23Meth [ 85 ]  =k32Meth[85]  * f 32 ( *Hz /Pa* ) ; 
k23Meth [ 87 ]  =k32Meth[87]  * f32 (*Hz/Pa*)  ; 
k32He [ 85 ]  =  1 . 1776  *  10s (*Hz/Pa*) ; 
k32He [ 87 ]  =  1 . 1776  *  10s (*Hz/Pa*) ; 
k23He [ 85 ]  =k32He[85]  * f32 (*Hz/Pa*) ; 
k23He [ 87 ]  =k32He[87]  * f32 (*Hz/Pa*) ; 
k32Alk [ 85 ]  =0(*Hz/Pa*); 
k32Alk [ 87 ]  =0(*Hz/Pa*); 
k23Alk [ 85 ]  =k32Alk[85]  * f 32 ( *Hz /Pa* ) ; 
k23Alk [ 87 ]  =k32Alk[87]  * f 32 ( *Hz /Pa* )  ; 

V.  Lineshape,  Emission,  and  Absorbtion  Cross  Sections 

A.  Transition  Lineshape 

(*Einstein  B  coeffs  based  on  the  A  coeff  and  the  frequncy  in  Hz*) 
c3  *  A21 [ iso] 

B21[v_,  iso_]  :=  - (*m3/j*s2  or  m/kg  which  are  the  same*); 

n2  *  ng  *  8  7r  *  h  *  v3 

B12[v_,  iso__]  :  =  g2  /  gl  *  B21  [  v,  iso]  (*m3/j*s2  or  m/kg*); 
c3  *  A31 [ iso] 

B31[v_,  iso__]  :=  - (*m3/J*s2  or  m/kg*); 

n2  *  ng  *  8  7r  *  h  *  v3 

B13[v_,  iso_]  :  =  g3  /  gl  *  B31  [  v,  iso]  (*m3/j*s2  or  m/kg*); 
rrad21  [iso_]  :=  (A21 [iso] ) _1 (*s*)  ; 
rrad31  [iso_]  :=  (A31 [iso] ) _1 (*s*)  ; 
rrad32  =  co(*s*)  ; 

(*Lorentzian  Linewidth  for  transitions*) 

1  1  I  tempMeth  1/2 

Avl21  [iso__]  :=  -  +  -&Methl  *  MMeth  *  -  + 

2  7r  rrad21  [iso]  (  temp 

( tempHe  1/2  tempAlk  1/2 

yHel  *  MHe  *  -  +  yAlkl  *  MAlk  *  -  (*Hz*)  ; 

(  temp  temp 

1  1  { tempMeth  1/2 

Avl31  [iso_]  :=  -  + ->Meth2  *  MMeth  *  -  + 

2  n  rrad31  [iso]  (  temp 

I  tempHe  1/2  tempAlk  1/2 

yHe2  *  MHe  *  -  +  yAlk2  *  MAlk  *  -  (*Hz*)  ; 

(  temp  temp 

Mass [85]  =  84 . 911789738  * amu (*Lewis  in  kg*); 

Mass [87]  =  86.909180527  *  amu  (*Lewis  in  kg*); 
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(*Hyperfine  corrections  to  frequency  for  87  and  85  rb  indexed 
by  F'  '  for  low  and  F'  for  ups.  Note  :  561d  is  F '  '  or  Sl/2  level*) 
vhyD21ow [ 2 ,  87]  =  2 . 563005979089114  * 109 (*Hz*) ; 
vhyD21ow [ 1 ,  87]  =- 4 . 271676631815196  * 109 (*Hz*) ; 
vhyD2up  [3,  87]  =  193 . 740846  *  106 (*Hz*)  ; 
vhyD2up  [2 ,  87]  =  72 . 911332  *  106 (*Hz*) ; 
vhyD2up  [  1 ,  87]  =- 229 . 851856  *  106 ( *Hz*) ; 
vhyD2up  [ 0 ,  87]  =- 302 . 073888  *  106 ( *Hz*) ; 
vhyDllow  [ 2 ,  87]  =  2 . 563005979089114  *  109 (*Hz*) ; 
vhyDllow  [  1 ,  87]  =- 4 . 271676631815196  *  109 (*Hz*) ; 
vhyDlup  [ 2 ,  87]  =  306 . 246  *  106 (*Hz*) ; 
vhyDlup  [  1 ,  87]  =  -  510 . 410  * 106 (*Hz*)  ; 

vhyD21ow  [3,  85]  =  1 . 26488516325  *  109 (*Hz*) ; 
vhyD21ow [ 2 ,  85]  =- 1 . 770843922835  * 109 (*Hz*) ; 
vhyD2up  [ 4 ,  85]  =  100 . 20544  *  106 (*Hz*) ; 
vhyD2up  [3,  85]  =  20 . 43551  *  106 (*Hz*) ; 
vhyD2up  [ 2 ,  85]  =- 83 . 83534  *  106 ( *Hz*) ; 
vhyD2up  [  1 ,  85]  =- 113 . 20884  *  106 (*Hz*)  ; 
vhyDllow  [3,  85]  =  1 . 26488516325  *  109 (*Hz*) ; 
vhyDllow [2,  85]  =- 1 . 770843922835  * 109 (*Hz*) ; 
vhyDlup  [3,  85]  =  150 . 65971  *  106 (*Hz*) ; 
vhyDlup  [2,  85]  =  -  210 . 4 92399  *  106 ( *Hz*) ; 
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(*Line  positions  for  hyperfine  lines  need  several  enteries  1,2,3, 

4  for  Sews  and  561ds  (These  are  usually  F'  and  F' '  in  literature  where  F  is  I  , 
the  nucler  spin,  +  J,  the  angular  momentum*) 

(*vhy21  [®ld  ,  Sew  ,  £>,  iso]  =  :  ;  *) 

(*vhy31  [®ld  ,  5fiew  ,  £>,  iso]  =  :  ;  *) 


vhy21  [561d_  ,  5fiew_  ,  ,  iso_] 

vhy31  [561d_  ,  5fiew_  ,  £>_  ,  iso_] 

A\d21  [561d_  ,  5fiew_  ,  ,  iso_] 

vhy21  [561d  ,  sfiew  ,  £>,  iso]  * 


=  v21  -  (vhyDllow  [561d  ,  iso]  )  +  vhyDlup  [5Rew  ,  iso]  (*Hz*)  ; 
=  v31  -  (vhyD21ow  [561d  ,  iso]  )  +  vhyD2up  [5Rew  ,  iso]  (*Hz*)  ; 


8  *  kb  *  tempcell  *  Log  [  2  ] 


Mass [iso]  *  c2 


A\d31  [561d_  ,  5Rew__  ,  ,  iso_]  :=  vhy31  [561d  ,  5fiew  ,  £>,  iso]  * 

(*Hz*)  ; 

a21[v_,  ®ld„  ,  Sew_  ,  B_  ,  iso_]  :=  (Log[2])2  * 


(*Hz*) ; 

8  *  kb  *  tempcell  *  Log [ 2 ] 


Mass [iso]  *  c2 


Avl21  [iso] 


(*Unitless*) ; 


Avd21  [551d  ,  Slew  ,  £>,  iso] 

1 

u21  [  v_  ,  !61d_  ,  5Rew_  ,  £>_  ,  iso_]  :  =  2  *  (Log  [2  ]  )  I  * 

|v  -  vhy21  [561d  ,  Sew  ,  £>,  iso]  (*  +  dMel  *MMeth*  ^  te”pMath  j  +6Hel  *MHe*  ^  t°mpHa  |  *)jy/ 

Avd21  [5Bld  ,  5fiew  ,  2),  iso]  (*Unitless*)  ; 

z21[v_,  !Sld_  ,  5Rew_  ,  ,  iso_]  :=  a21[v,  ®ld  ,  5fiew  ,  D ,  iso]  +i*u21[v,  ®ld  ,  Shew  ,  D,  iso] 

(*Unitless*) ; 


a31[v_,  !Sld_  ,  Sew_  ,  ,  iso_J  :=  (Log  [2])  2  * 


Avl31  [iso] 


(*Unitless*) ; 


Avd31  [561d  ,  Slew  ,  £>,  iso] 

1 

u31  [  v_ ,  ®ld_  ,  Sew_  ,  ,  iso_]  :  =  2  *  (Log  [2  ]  )  2  * 

(/  tempMeth  \  ^ ^ ^  /  tempHe  \  \  / 

v  -  vhy31  [561d  ,  Sew  ,  D,  iso]  (*  +  5Me2  *MMeth*  | — - J  +<5He2  *MHe*  |  — - j  *)  )  / 

Avd31  [561d  ,  5fiew  ,  D,  iso]  (*Unitless*)  ; 

z31[v_,  !Sld_  ,  5Rew_  ,  ,  iso_]  :=  a31[v,  ®ld  ,  5Rew  ,  D ,  iso]  +i*u31[v,  ®ld  ,  5fiew  ,  B,  iso] 

(*Unitless*) ; 

(*The  lineshape  of  each  transition.  These  functions  represent 
sets  of  functions  for  each  transition  indexed  by  there  hyperfine 
compenents  F'  and  F'  '  and  D  the  manifold  or  5fiew  and  561d  and  D*) 


g21[v_,  !Sld_  ,  Sew_  ,  1,  iso_]  :  =  NIntegrate 


(  Mass [iso] 


^  2  n  *  kb  *  temp  J 

Avl21  [iso]  /  (2  7r  *  (  (v  -  vhy21  [5Bld  ,  Sew  ,  1,  iso]  -  vhy21  [5Bld  ,  Sew  ,  1,  iso]  *  vz  /  c)  2  + 
( Avl21  [iso]  /  2)  2)  )  ,  {vz  ,  -106,  106}]/  1012  (*  in  s*)  ; 


g31[v_,  !Sld__  ,  5Rew_  ,  2,  iso_]  :=  NIntegrate 


Mass [iso] 


^  2  7r  *  kb  *  temp  , 

Avl31  [iso]  /  (2  7T  *  (  (v  -  vhy31  [5Bld  ,  3iew  ,  2,  iso]  -  vhy31  [561d  ,  Shew  ,  2,  iso]  *  vz  /  c)  2  + 
( Avl31  [iso]  /  2)  2)  )  ,  {wz,  -1010,  1010}]/  1012  (*in  s*)  ; 
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g21approx  [  v_  ,  !Bld_  ,  3iew_ 

rs  |^(e(z21[v,^ld  ,®ew  ,  1,  iso]  )  2 

g31approx  [  v_  ,  ®ld_  ,  31ew_ 


,  1,  iso_]  := 


4  Log [ 2 ] 


Erfc[z21[v,  561d  ,  5Rew 
f 4  Log [ 2 ] 

,  2,  iso_]  := 


Re[e<z31[v':51d  '*ew  ■2-is°])2  Erfc[z31  [v,  561d  ,  5Rew  , 


1/2 


1, 

1/2 


2, 


1 

Avd21  [561d  ,  SRew  ,  1,  iso] 
iso] ] ]  ( * s * )  ; 

1 

Avd31  [®ld  ,  fiew  ,  2,  iso] 
iso] ] ]  ( * s * ) ; 


gapproxfix21  [  v_ ,  ®ld_  ,  8ew_  ,  1,  iso_]  :  = 

If  [g21approx  [  v,  5Bld  ,  3iew  ,  1,  iso]  =  =  =  Underflow  []  , 

0,  g21approx[v,  ®ld  ,  8ew  ,  1,  iso]]  (*s*) ; 
gapproxfix31  [  v_ ,  ®ld_  ,  JRew_  ,  2,  iso_]  :=  If  [g31approx  [  v,  ®ld  ,  8ew  ,  2,  iso]  =  =  = 
Underflow]],  0,  g31approx[v,  5Bld  ,  Slew  ,  2,  iso]]  (*s*)  ; 

(*Plot [ {gapproxf ix21 [v, 1,1, 1,87] , gapproxfix21 [v, 1,2, 1,87], gapproxfix21 [v, 2, 2, 1,87], 
gappr oxf ix2 1 [v,2,l,l,87]},{v, vhy21 [2, 2, 1, 87] +2*109, vhy21 [ 1 , 1 , 1 , 87 ] -2*109 } , 
PlotRange->{  { vhy21  [1, 2, 1, 87]  -109,  vhy21  [2, 1, 1, 87]  +109}  , 

{ 0 , gapproxf ix21 [vhy21 [1,2,1,87] , 2 , 2 , 1 , 87 ] *2} } ] 

Plot [ {gapproxfix31 [v, 1,0, 2, 87] , gapproxf ix31 [v, 1,1, 2, 87], gapproxf ix31 [v, 1,2, 2, 87] , 
gapproxfix31 [v, 2,2,2,87] , gapproxfix31 [v, 2,1,2, 87] , gapproxfix31 [v, 2,3,2, 87] } , 

{ v , vhy 3 1  [2,2,2, 87] -2*109, vhy31  [1,  1, 2, 87] +2*109}, 

PlotRange->{ { vhy31  [2, 2, 2, 87] +2*109, vhy31  [ 1 , 1 , 2 , 87 ] -2*109 } , 

{ 0 , gapproxf ix31 [vhy31 [2,2,2,87] , 2 , 2 , 2 , 87 ] *2} } ] *) 
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B.  Absorbtion  Cross  Section 
i.  Hyperfine  Structure 

(*.S[561d  ,  5fiew  ,D,iso]  is  the  relative  line  strength  All  are  unitless*) 
1 

■S[l,  1,  1,  87]  =  —  (*Unitless*)  ; 

6 

5 

■S[l,  2,  1,  87]  = 

6 
1 

■S  [2 ,  2,  1,  87]  = 

2 

1 

■S[2,  1,  1,  87]  = 

2 

1 

■S[l,  0,  2,  87]  = 

6 

5 

•S[l,  1,  2,  87]  =  — ; 

12 

5 

■S[l,  2,  2,  87]  =  — ; 

12 

■S[l,  3,  2,  87]  =  0; 

-S  [ 2 ,  0,  2,  87]  =  0; 

1 

■S  [2 ,  1,  2,  87]  =  — ; 

20 

1 

■S  [2 ,  2,  2,  87]  = 

4 
7 

■S  [ 2 ,  3,  2,  87]  =  — ; 

10 

2 

■S  [2 ,  2,  1,  85]  = 

9 

7 

■S  [  2 ,  3,  1,  85]  = 

9 

5 

■S  [  3 ,  2,  1,  85]  = 

9 

4 

■S  [  3 ,  3,  1,  85]  = 

9 
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S[ 2,  1,  2, 
5[2,  2,  2, 

5[2,  3,  2, 

5[2,  4,  2, 
5[3,  1,  2, 

5[3,  2,  2, 
S[3,  3,  2, 
5[3,  4,  2, 


3 

85]  =  — ; 
10 

7 

85]  =  — ; 
18 

14 

85]  =  — ; 
45 

85]  =  0; 
85]  =  0; 

5 

85]  =  — ; 
63 

5 

85]  =  — ; 
18 

9 

85]  =  — ; 
14 


Sf-,f1 

87Rb 

85Rb 

F 

1 

2 

F 

2 

3 

Dl  2S«2  -»  2P|I2 

1 

1/6 

1/2 

2 

2/9 

5/9 

2 

5/6 

1/2 

3 

7/9 

4/9 

D2  2S((2  — >  2Pai2 

0 

1/6 

... 

1 

3/10 

... 

1 

5/12 

1/20 

2 

7/18 

5/63 

2 

5/12 

1/4 

3 

14/45 

5/18 

3 

... 

7/10 

4 

... 

9/14 

(*hyperfine  line  strengths  for  structure 

need  them  for  each  forbidden  transitions  S  =  0*) 

Plot  [  {gapproxf  ix21  [  v,  1,  1,  1,  87]  *5[1,  1,  1,  87]  +  gapproxfix21  [v,  1,  2,  1,  87]  *S[1,  2,  1,  87]  + 

gapproxfix21  [v,  2,  2,  1,  87]  *S[2,  2,  1,  87]  +  gapproxfix21  [ v,  2,  1,  1,  87]  *d[2,  1,  1,  87]}, 

{v,  vhy2 1  [ 2 ,  2,  1,  87]  +PumpFWHM,  vhy21[l,  1,  1,  87]  -PumpFWHM}, 

PlotRange  ->  {{vhy21[2,  2,  1,  87]  +  PumpFWHM,  vhy21[l,  1,  1,  87]  -PumpFWHM}, 

{0,  gapproxf ix21  [vhy21  [2,  2,  1,  87],  2,  2,  1,  87]  *1}}, 

PlotLabel  -»  "g2i(v)  vs.  Freq .  For  87Rb",  AxesLabel  ->  {Frequency  [Hz],  g[v]  [unitless]}] 

Plot  [  {gapproxf  ix21  [v,  2,  2,  1,  85]  *d[2,  2,  1,  85]  +  gapproxf  ix21  [v,  2,  3,  1,  85]  *S[2,  3,  1,  85]  + 

gapproxf ix21  [v,  3,  3,  1,  85]  *S[3,  3,  1,  85]  +  gapproxfix21  [ v,  3,  2,  1,  85]  *d[3,  2,  1,  85]}, 

{v,  vhy2 1  [ 3 ,  3,  1,  85]  +  PumpFWHM,  vhy21[2,  2,  1,  85]  -  PumpFWHM}, 

PlotRange  ->  {{vhy21[3,  3,  1,  85]  +PumpFWHM,  vhy21[2,  2,  1,  85]  -PumpFWHM}, 

{0,  gapproxf ix21  [vhy21  [2,  2,  1,  85],  2,  2,  1,  85]  *1}}, 

PlotLabel  -»  "g2i(v)  vs.  Freq.  For  85Rb",  AxesLabel  ->  {Frequency  [Hz],  g[v]  [unitless]}] 

Plot  [  {gapproxf  ix31  [  v,  1,  0,  2,  87]  *d[l,  0,  2,  87]  +  gapproxfix31  [v,  1,  1,  2,  87]  *S[1,  1,  2,  87]  + 

gapproxfix31  [v,  1,  2,  2,  87]  *S[1,  2,  2,  87]  +  gapproxf ix31  [v,  2,  2,  2,  87]  *,S[2,  2,  2,  87]  + 
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gapproxfix31  [v,  2,  1,  2,  87]  *S[2,  1,  2,  87]  +  gapproxfix31  [ v,  2,  3,  2,  87]  *£[2,  3,  2,  87]}, 
{v,  vhy 3 1  [ 2 ,  2,  2,  87]  +PumpFWHM,  vhy31[l,  1,  2,  87]  -PumpFWHM}, 

PlotRange  ->  {{vhy31[2,  2,  2,  87]  +  PumpFWHM,  vhy31[l,  1,  2,  87]  -  PumpFWHM}, 

{0,  gapproxfix31  [vhy31  [2,  2,  2,  87],  2,  2,  2,  87]  *1}}, 

PlotLabel  -»  "g3i(v)  vs.  Freq .  For  87Rb",  AxesLabel  ->  {Frequency  [Hz],  g[v]  [unitless]}] 

Plot  [  {gapproxf  ix31  [  v,  2,  1,  2,  85]  *5[2,  1,  2,  85]  +  gapproxfix31  [v,  2,  2,  2,  85]  *S[2,  2,  2,  85] 

gapproxfix31  [v,  2,  3,  2,  85]  *S[2,  3,  2,  85]  +  gapproxf ix31  [v,  3,  3,  2,  85]  *,S[3,  3,  2,  85]  + 

gapproxfix31  [v,  3,  2,  2,  85]  *S[3,  2,  2,  85]  +  gapproxf ix31  [v,  3,  4,  2,  85]  *,S[3,  4,  2,  85]}, 

{v,  vhy 3 1  [ 3 ,  3,  2,  85]  +PumpFWHM,  vhy31[2,  2,  2,  85]  -PumpFWHM}, 

PlotRange  ->  {{vhy31[3,  3,  2,  85]  +  PumpFWHM,  vhy31[2,  2,  2,  85]  -  PumpFWHM}, 

{0,  gapproxf ix31  [v31  ,  3,  3,  2,  85]  *15}}, 

PlotLabel  -»  "g3i(v)  vs.  Freq.  For  85Rb",  AxesLabel  ->  {Frequency  [Hz],  g[v]  [unitless]}] 
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ii.  Statistical  isotope  distribution 


fiso[85]  =  0 . 7217 ( *Unitless *) ; 
fiso[87]  =  0 . 2783  (*Unitless *)  ; 

SRewmin  [1,  87]  =  1 (*Unitless*) ; 
SRewmax  [1,  87]  =  2 ( *Unitless *) ; 
SBldmin  [1,  87]  =  1 (*Unitless*) ; 
SBldmax  [1,  87]  =  2 ( *Unitless *) ; 

SRewmin  [2,  87]  =  0 (*Unitless*) ; 
SRewmax  [2,  87]  =  3 ( *Unitless *) ; 
SBldmin  [2,  87]  =  1 (*Unitless*) ; 
SBldmax  [2,  87]  =  2 ( *Unitless *)  ; 

SRewmin  [1,  85]  =  2 (*Unitless*) ; 
SRewmax  [1,  85]  =  3 ( *Unitless *) ; 
SBldmin  [1,  85]  =  2 (*Unitless*) ; 
SBldmax  [1,  85]  =  3 ( *Unitless *) ; 

SRewmin  [2,  85]  =  1 (*Unitless*) ; 
SRewmax  [2,  85]  =  4 ( *Unitless *) ; 
SBldmin  [2,  85]  =  2 (*Unitless*) ; 
SBldmax  [2,  85]  =  3 ( *Unitless *) ; 

isomin  =  85 (*Unitless*) ; 
isomax  =  87 (*Unitless*) ; 
isostep  =  2 (*Unitless*) ; 


EngD2  [561d_  ,  iso_]  :  =  h*  (v31  -  vhyD21ow  [SBld  ,  iso]  ) 

(*The  energy  of  a  given  F''  state  in  J*) ; 

EngDl  [561d_  ,  iso_]  :  =  h  *  (v21  -  vhyDllow  [SBld  ,  iso]  )  (*J*)  ; 

®ldmax  [1,85]  -EngDl  [sold  ,  851 

ftotDl[85]  =  ^  (2  SBld  +  1)  *  e  kb"t"r*'  (*Unitless*)  ; 

5Bld  =5Bldmin  [1,85] 

[©ldmax  [1,87]  -EngDl  [sold  ,  87] 

ftotDl[87]  =  ^  (2  SBld  +  1)  *  e  kb"t"rf  (*Unitless*)  ; 

5Bld  =5Bldmin  [1,87] 

®ldmax  [2,85]  -EngD2  [SBld  ,  85] 

ftotD2[85]  :=  ^  (2  SBld  +  1)  *  e  k“’*t"r  (*Unitless*)  ; 

561d  =561dmin  [2,85] 

^ldmax  [2,87]  -EngD2  [SBld  ,  87] 

ftotD2[87]  :=  ^  (2  SBld  +  1)  *  e  kb*t“I1P  (*Unitless*)  ; 

5Bld  =561dmin  [2,87] 

-EngDl  [SBld  ,  iso] 

(2  SBld  +  1)  *  <E~  kb* temp 

ff21[561d_  ,  iso_]  :=  - 

ftotDl [ iso] 

(*stat  dis  among  T  states  will  need  several  enteries  for  each (*Unitless *) *) 

-EngD2  [SBld  ,  iso] 

(2  SBld  +  1)  *  e  kb* temp 

ff31[5Bld_  ,  iso_]  :=  - 

ftotD2 [ iso] 

(*stat  dis  among  T  states  will  need  several  enteries  for  each (*Unitless*) *) 
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iii.  Cross  Sections  (in  m2) 


al2abs [v_]  = 

1  561dmax  [D,  iso]  Jiewmax  [D,iso] 


sum[Z  z 

D=  1  561d  =5Bldmin  [D,  iso]  Slew  =2ievnnin  [D,  iso] 


z 


g2 

U1 


t  V2  *  8  7T 


A21[iso]  *5[®ld  ,  Slew  ,  £>,  iso]  * 


gapproxfix21  [v,  561d  ,  5Rew  ,  £>,  iso]  *  ff21  [5Bld  ,  iso]  *fiso[iso] 
{iso,  isomin,  isomax,  isostepjj  (*m2*) ; 


2  561dmax  [D,  iso]  Jiewmax  [D,iso] 

ol3abs  [v_  ]  =  Sum[Z  Z  Z 

Z)=2  561d  =5Bldmin  [D,  iso]  Jiew  =2iewmin  [Z),iso] 


g3 

U1 


\  V2  *  8  7T 


A31[iso]  *d[®ld  ,  3 lew  , 


iso]  *  gapproxfix31  [v,  ®ld  ,  fiew  ,  £>,  iso]  *  ff31  [®ld  ,  iso]  *fiso[iso] 
,  {iso,  isomin,  isomax,  isostepjj  (*m2*) ; 
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C.  Emission  Cross  Sections 


gl 

o21ems  [ v _ ]  =  —  al2abs  [v]  (*m2*)  ; 

g2 


gl  5 

cr31ems  [v_]  =  —  ol3abs  [v]  (*m  *)  ; 

g3 


(*Do [print [o21ems  [vhy21  [2,2,1,87]  -106  +  i*105]  ]  ,  {i,  100}]  ;  *) 
(*trapend  =  2*101:1 ; 
trapstep=109 ; 

trapstart  =vhy21  [2,2,1,87]-109; 


trapend 

trapend/trapstep 


*Sum [o21ems  [vhy21  [2,2,1,87]  -101:L  +  l]  ,  {1,0,  trapend,  trapstep}  ]  *) 


(*Plot [{o21ems [v_]  }  , {v,  vhy21  [2,2,1,87] -109, vhy21 [2, 2, 1, 87] +109}  , 

PlotRange-»{ { vhy21  [2, 2, 1, 87]  -109,  vhy21  [2 , 2 , 1 , 87 ]  +109}  ,  { 0 , 2*10-7 }  }  ]  *) 

(*NIntegrate [o21ems [ v_] , { v, vhy21 [2,2,1,87] , vhy21 [1,1,1,87]}, Method-»"TrapezoidalRule " ] *) 
gapproxfixtotal21 [v_ ]  = 

1  561dmax  [D,  iso]  Jiewmax  [D,  iso] 

Sum[ ^  ^  gapproxfix21  [v,  561d  ,  Bew  ,  D,  iso]  * 

Z)=l  561d  =561dmin  [D,  iso]  Jiew  =3iewmin  [2),iso] 

>S[®ld  ,  3iew  ,  D,  iso]  *fiso[iso]  *  ff21  [®ld  ,  iso],  {iso,  isomin,  isomax,  isostep}j; 

2  SBldmax  [Z),  iso]  Jiewmax  [Z),iso] 

gapproxfixtotal31  [v_]  =  Sum[^  ^  ^  gapproxfix31  [ 

D=2  561d  =5Bldmin  [Z),iso]  Slew  =®ewmin  [Z),iso] 

v,  ®ld  ,  3iew  ,  £>,  iso]  *d[®ld  ,  Slew  ,  £),  iso]  * 
fiso[iso]  ff31  [®ld  ,  iso],  {iso,  isomin,  isomax,  isostep}j; 

Plot  [  {gapproxf  ixtotal31  [  v31  +  v]  ,  glpO  [  v31  +  v]  }  ,  {v,  -2*PumpFWHM,  2*PumpFWHM}, 

PlotRange  {  {  -  PumpFWHM  *  2 ,  PumpFWHM*2},  {0,  gapproxf  ixtotal31  [  vhy31  [  2 ,  2,  2,  87]]  *3}}, 
PlotLabel  ->  Style  [  "g31  (v)  (Blue)  and  gpump(vp)  (Red)  for  Rb  vs.  Freq .  " ,  24], 

LabelStyle  -> Directive [Large] ,  AxesLabel  ->  {Frequency  [Hz],  g[v]  [unitless]}, 

LabelStyle  ->  Directive [Large] ,  AxesOrigin  ->  {-PumpFWHM  *  2,  0}] 

Plot [o31ems  [v31  +  v]  ,  {v,  -2  *  PumpFWHM,  2  *  PumpFWHM} , 

PlotRange  ->  {{- 2  *  PumpFWHM ,  2*PumpFWHM},  {0,  o31ems  [  vhy31  [2 ,  2,  2,  87]]  *3}}, 

PlotLabel  ->  Style [ "a3iems [v]  vs.  Freq.",  24],  LabelStyle  -»  Directive [Large] , 

AxesLabel  ->  {  Frequency  [  Hz]  ,  o31>ms  [v]  [m2]  ] ,  LabelStyle  -»  Directive  [Large]  , 

AxesOrigin  ->  {  -  2  * PumpFWHM,  0},  PlotPoints  ->  400] 


Plot [gapproxfixtotal21  [  v21  +  v] ,  {v,  -  1  /  5  *  PumpFWHM,  1  / 4  *  PumpFWHM} , 

PlotLabel  ->  Style  [”g2i  (v)  For  Rb  vs.  Freq.",  24], 

AxesLabel  -»  {Frequency  [Hz],  g[v]  [unitless]},  LabelStyle  -» Directive [Large] , 

AxesOrigin  ->  {-PumpFWHM  *1/5,  0}  ,  PlotRange  -* 

{  {-PumpFWHM  *  1  /  5,  PumpFWHM  *  1  /  4 }  ,  {0,  gapproxf ixtotal21  [ vhy21  [  1 ,  1,  1,  87]]  *10}}] 
Plot  [o21ems  [v21  +  v]  ,  {v,  -2*  PumpFWHM,  +  2  *  PumpFWHM }  , 

PlotRange  -»  {{ -PumpFWHM  *  2 ,  PumpFWHM  *2},  {0,  o21ems  [ vhy21  [  1 ,  1,  1,  87]]  *6}}, 

PlotLabel  ->  Style  [  "o2iems  (v)  vs.  Freq.",  24],  AxesLabel  ->  {  Frequency  [  Hz  ]  ,  u21m-  [  v]  [m2  ]  }  , 
LabelStyle  ->  Directive [Large] ,  AxesOrigin  -»  {-PumpFWHM  *  2,  0},  PlotPoints  ->  400] 
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D.  Cavity-Mode  Spacing 

c 

AvEsr  =  - ; 

2  *  n  *  distance 

vmaxgain  =  v  /  .  Last  [FindMaximum  [o21ems  [v]  ,  {v,  v21  -  1010,  v21  +  1010}  ]  ]  ; 
qmaxgain  =  Round [vmaxgain  /  AvEsr  ]  ; 
vlaser  =  qmaxgain  *  AvEsr ; 

r  vlaser  o21ems  [v21  +  v] 

Plot - ,  {v,  -  2  *  PumpFWHM ,  2  *  PumpFWHM}  , 

Pumplinecenter  o31ems  [Pumplinecenter  +  v] 

PlotLabel  ->  Style ["rj  The  effciency  factor  vs.  Frequency  in  Hz",  24], 

AxesLabel  ->  [Freq[Hz]  ,  77 }  ,  LabelStyle  ->  Directive  [Large]  1 
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VII.  Rate  Equations  Analysis,  ODE  Solutions,  and  System  Outputs 


) 
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A.  Species  Setups 

(*If  you  wish  to  add  more  species  you  must  do  so  in  II.  C.  , 

I.  D.,  and  in  V.  A.  Also  rember  the  species  numer  (1,2,3...) 
must  be  the  same  for  all  coeffs  and  pressures .  [2]  must  always 

correspond  to  Heane  for  the  species,  but  other  sums  do  exist*) 
k23[l,  iso_]  :  =  k23Meth [iso]  (*Hz*) ; 
k23[2,  iso_]  :  =  k23He [iso]  (*Hz*)  ; 
k23[3,  iso_]  :  =  k23Alk [iso]  (*Hz*) ; 
k32[l,  iso_]  :  =  k32Meth [iso]  (*Hz*) ; 
k32[2,  iso_]  :  =  k32He [iso]  (*Hz*) ; 
k32[3,  iso_]  :  =  k32Alk [iso]  (*Hz*) ; 

M[l]  =  MMeth;  (*Partial  Pressures  in  Pa*) 

M  [2  ]  =  MHe ( *Pa* ) ; 

M [ 3  ]  =  MAlk (*Pa*) ; 
speciesmax  =  3 (*Unitless*) ; 
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B.  Rate  Equations  and  Populations 


M  =  TotalAlkaliConcentration  (*in  m  3*)  ; 

II [z_]  :  =  Iplu [z]  + Iminu [z]  (*w/m2*)  ; 

(*Reqnl[At_  ,J2_  ,12_  ]  :  =  *A31+/2  *A21  + 

12  *k31+j2  *k21+ (-B13  [v]  */I  *Ip+B31  [  v]  */3  *Ip+B21  [  v]  *11*  -^-/I  ))dlv; 

Reqn2  [/!_  ,  J2_  ,  12_  ]  :  =  JC^Zspacies=iX-^  *k32Meth*MMeth-/2  *k23Meth*MMeth- 
12  *A21  -J2  *k21  -B21  [v]  *11*  )dlv; 

Reqn3  [/!_  ,  12_  ,  12_  ]  :  =  JC^Zspacias=iX-/3  *k32Meth*MMeth+J2  *k23Meth*MMeth- 
!2  *A31-/3  *k31  +B13  [  v]  *Ip*/I  -B31  [v]  *12  *Ip  dlv; 

Reqn4[/I_  ,  12__  ,  12_  ]  :  =/I  +J2  +/S  ; 

Solve  [  {Reqn2  [12.  ,12  ,12  ],Reqn3[/I  ,  12  ,  12  ],Reqn4[/I  ,12  ,12  ]}=={0,0,/t£  },{/!  ,12  ,12  }]; 


- ;  *) 


J2  +/S 

Pop3[z_,  v_,  vpump_]  = 


Sum  Sum  fiso[iso] 


speciesmax 


JV. 


■  k21  -  A21 [iso]  -  k23 [ species ,  iso]  M[species] 
g2Il[z]  o21ems  [ vlaser  ] 


II [z]  o21ems [vlaser ] 


gl h  vlaser 


h  vlaser 

g2  II [z]  o21ems  [vlaser]  g3  Ip[z]  a31ems  [vpump] 


gl h  vlaser 


g2  M  II [z]  o21ems  [vlaser ] 


k21  +  A21  [iso]  + 


gl h  vpump 
Il[z]  o21ems [vlaser ] 


gl h  vlaser 

g2  II [z]  o21ems  [vlaser]  g3  Ip[z]  a31ems  [vpump] 


gl h  vlaser 

k32 [species,  iso]  M[species] 


gl h  vpump 
g2  II [z]  o21ems [vlaser] 


/ 


h  vlaser 


k21  +  A21  [iso]  + 


gl h  vlaser 

II  [z]  o21ems  [vlaser]  g2  Il[z]  o21ems  [vlaser]  g3  Ip[z]  a31ems  [vpump] 


h  vlaser  gl h  vlaser 

-k21  -  A21 [iso]  -  k23 [ species ,  iso]  M[species] 
g2  II [z]  o21ems  [vlaser] 


gl h  vpump 
II  [z]  o21ems  [vlaser] 


gl h  vlaser 


k31  +  A31  [iso]  + 


h  vlaser 
g2  II [z]  o21ems [vlaser] 


gl h  vlaser 
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Ip[z]  o31ems  [vpump]  g3  Ip[z]  o31ems  [vpump] 


h  vpump  gl h vpump 

{species,  1,  speciesmax,  1}]/  {iso,  isomin,  isomax,  isostep }  j  (*m~3*)  ; 

1 


Pop2[z_,  v_  ,  vpump_ ]  = 


Sum  Sum  fiso[iso] 


speciesmax 


(g2  M  II [z]  o21ems  [  vlaser  ]  )  /  gl hvlaser 


/(? 


- k21  -  A21 [iso]  - 


k23[species,  iso]  M[species]  - 

k32[ species,  iso]  M[ species]  - 


II  [z]  o21ems  [vlaser  ]  g2Il[z]  o21ems  [vlaser  ] 


h  vlaser 
g2Il[z]  o21ems  [vlaser  ] 


gl h  vlaser 


M 


gl h  vlaser 

-k21  -  A21 [iso]  -  k23 [ species ,  iso]  M[species]  - 
g2Il[z]  o21ems  [vlaser  ] 


II [z]  o21ems  [vlaser] 


gl h  vlaser 


h  vlaser 

g2  II [z]  o21ems  [vlaser]  g3  Ip[z]  o31ems  [vpump] 


gl h  vlaser 


g2  M  II [z]  o21ems  [vlaser ] 


k21  +  A21  [iso]  + 


gl  h  vpump 
II  [z]  o21ems  [vlaser] 


gl  h  vlaser 

g2  II  [z]  o21ems  [vlaser]  g3  Ip[z]  o31ems  [vpump] 


h  vlaser 


gl h  vlaser 


gl  h  vpump 


/ 


-k21  -  A21 [iso]  -  k23 [ species ,  iso]  M[species]  - 
g2  II [z]  o21ems [vlaser] 


II [z]  o21ems  [vlaser] 


h  vlaser 


gl h  vlaser 


k32[ species,  iso]  M[ species]  - 


g2  II [z]  o21ems [vlaser] 


k21  +  A21  [iso]  + 


gl  h  vlaser 

Il[z]  o21ems  [vlaser]  g2  Il[z]  o21ems  [vlaser]  g3  Ip[z]  o31ems  [vpump] 


h  vlaser  gl h  vlaser 

-k21  -  A21 [iso]  -  k23 [species ,  iso]  M[species] 
g2Il[z]  o21ems  [vlaser  ] 


gl  h  vpump 
Il[z]  o21ems [vlaser ] 


k31  +  A31  [iso]  + 


h  vlaser 
g2  II [z]  o21ems [vlaser] 


gl h  vlaser 

Ip[z]  o31ems  [vpump]  g3  Ip[z]  o31ems  [vpump] 


gl  h  vlaser 


h  vpump  gl  h  vpump 

{species,  1,  speciesmax,  l}j,  {iso,  isomin,  isomax, 
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isostep} J  (*m  3*)  ; 

Popl[z_,  v_ ,  vpump_  ]  =  Sumj^Sum  Jfiso  [iso] 

( M  +  (g2  M  Il[z]  o21ems  [  vlaser  ]  )  /  I  gl  h  vlaser 


speciesmax 


'/(' 


-  k21  -  A21 [iso] 


k23[species,  iso]  M[species]  - 


II  [z]  o21ems  [vlaser  ]  g2Il[z]  o21ems  [vlaser  ] 


M  | -  k21  -  A21 [iso]  - k23 [ species ,  iso]  M [species] 
g2Il[z]  o21ems  [vlaser  ] 


h  vlaser  gl h  vlaser 

II [z]  o21ems [vlaser ] 


gl h  vlaser 


h  vlaser 

g2  II [z]  o21ems  [vlaser]  g3  Ip[z]  a31ems  [vpump] 


gl h  vlaser 


g2  M  II [z]  o21ems  [vlaser ] 


k21  +  A21  [iso]  + 


gl h  vpump 
Il[z]  o21ems [vlaser ] 


gl h  vlaser 

g2  II [z]  o21ems  [vlaser]  g3  Ip[z]  o31ems  [vpump] 


gl h  vlaser 

k32[ species,  iso]  M[ species] 


gl h  vpump 
g2Il[z]  o21ems  [vlaser  ] 


/ 


h  vlaser 


k21  +  A21  [iso]  + 


gl h  vlaser 

Il[z]  o21ems  [vlaser  ]  g2Il[z]  o21ems  [vlaser  ]  g3Ip[z]  o31ems  [vpump  ] 


h  vlaser  gl h  vlaser 

- k21  -  A21 [ iso]  -  k23 [species ,  iso]  M[species]  - 
g2Il[z]  o21ems  [vlaser  ] 


k31  +  A31  [iso]  + 


gl  h  vpump 
Il[z]  o21ems  [vlaser ] 
h  vlaser 
g2Il[z]  o21ems  [vlaser  ] 


gl h  vlaser 

Ip[z]  a31ems  [vpump]  g3  Ip[z]  o31ems  [vpump] 


gl h  vlaser 


h  vpump 

k32[ species,  iso]  M[ species] 


gl h  vpump 
g2Il[z]  o21ems  [vlaser  ] 


M 


gl  h  vlaser 

-k21  -  A21 [iso]  -  k23 [ species ,  iso]  M[species] 
g2Il[z]  o21ems  [vlaser  ] 


II [z]  o21ems  [vlaser] 


gl h  vlaser 


h  vlaser 

g2  II [z]  o21ems  [vlaser]  g3  Ip[z]  o31ems  [vpump] 


gl h  vlaser 


gl  h  vlaser 


g2  M  II [z]  o21ems  [vlaser ] 


k21  +  A21  [iso]  + 


gl  h  vpump 
II  [z]  o21ems  [vlaser] 


h  vlaser 
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g2  II  [z]  o21ems  [vlaser]  g3  Ip[z]  o31ems  [vpump] 
-  +  - 


gl h  vlaser 


gl  h  vpump 


)/ 


-k21  -  A21 [iso]  -  k23 [ species ,  iso]  M[species]  - 
g2  II [z]  o21ems [vlaser] 


II [z]  o21ems  [vlaser] 


h  vlaser 


gl h  vlaser 


k32 [ species,  iso]  M[ species] 


g2  II [z]  o21ems [vlaser] 


k21  +  A21  [iso]  + 


gl  h  vlaser 

Il[z]  o21ems  [vlaser]  g2  Il[z]  o21ems  [vlaser]  g3  Ip[z]  o31ems  [vpump] 


h  vlaser  gl h  vlaser 

-k21  -  A21 [iso]  -  k23 [species ,  iso]  M[species] 
g2Il[z]  o21ems  [vlaser  ] 


gl h  vlaser 


k31  +  A31  [iso]  + 


gl h  vpump 
Il[z]  o21ems [vlaser ] 
h  vlaser 
g2  II [z]  o21ems [vlaser] 


gl  h  vlaser 


Ip[z]  a31ems  [vpump]  g3  Ip[z]  o31ems  [vpump] 
-  +  - 


h  vpump  gl  h  vpump 

{species,  1,  speciesmax,  1 }  j  ,  {iso,  isomin,  isomax, 
isostep}  j  (*m“3*)  ; 


C.  f(stat  dis)  and  Differential  Equations  in  z  for  l+,  L,  and  L 


I sat [vpump_ ] 

fpop[z_,  v_] 


h  *  vpump 

-  i 

o31ems  [  vpump  ] 
Pop2 [z,  v] 


*  fiso [85]  + 


Pop2[z,  v]  +Pop3[z,  v] 


trad31  [85] 

(*Unitless*) ; 


rrad31 [87] 


*  fiso [87] 


(*Wm  2*)  ; 


a  =  -  *  Log  -  ;  ( *Loss  Coef f  *) 

2  *  celllength  rl  *  r2  *  tl  A  4  ^ 

Deqnl [z_,  v_  ,  vpump_  ]  :  = 

o31ems  [vpump  ]  ((Pop3[z,  v,  vpump])  -2*Popl[z,  v,  vpump]  )  *  Ip  [z]  (*m  1*)  ; 

Deqn2  [z_,  v_,  vpump_  ]  :=  (o21ems  [v]  ((Pop2[z,  v,  vpump])  -Popl[z,  v,  vpump]))  *Iplu[z] 
(*m-1*)  ; 

Deqn3  [z_,  v_,  vpump_  ]  :=  -  (o21ems  [v]  ((Pop2[z,  v,  vpump])  -Popl[z,  v,  vpump]))  *  Iminu  [  z  ] 
(♦nr1*)  ; 

(*lp[0] =20000000 ; 

Iminu [0]=1750000; 

Iplu [ 0 ] =1750000; 
vpump  =v31 ; 

vlaser  =vhy21 [2,2,1,87]  ; 

gapproxfixtotal211oop =gapproxf ixtotal21 [vlaser] ; 
glp01oop=glp0 [vpump ] ; 

Popl [ 0 , vhy21 [2,2,1,87] , vpump ] ; 

Pop2 [ 0 , vhy2 1 [2,2,1,87] , vpump ] ; 

Pop3 [ 0 , vhy21  [2,2,1,87] , vpump ]  ; 

Pop3 [ 0 , vhy21  [2,2,1,87],  vpump ] -2*Popl [ 0 , vhy21  [2,2,1,87], vpump ] ; 

Pop2 [ 0 , vhy21  [2,2,1,87] ,  vpump ] -Popl [ 0 , vhy21  [2,2,1,87], vpump ] 

Popl [ 0 , vhy2 1 [2,2,1,87] , vpump ] + 

Pop2 [0, vhy21  [2,2,1,87]  ,  vpump ] +Pop3 [ 0 , vhy21  [2,2,1,87], vpump ]  ; *) 
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D.  NDSolve  Routine  for  l+,  /_,  and  lp 

(*Ip  [  znot  ]  ==IpO*gIpO  [  v31  ]  ,  Iplu  [  znot  ]  ==Iminu  [  znot  ]  , 

refl*Iplu [ z final ]  ==Iminu [ z final] , Iplu [z] *Iminu [ z] ==Constant *) 

tl  =  cell-transmission  (*Unitless*)  ; 

rl  =  hrreflectivty ( *Unitless *) ; 

r2  =  ocreflectivty ( *Unitless *) ; 

znot  =  0.0 (*m*)  ; 

zfinal  =  znot  +  celllength (*m*) ; 

vrange  =  3  *  PumpFWHM ; 
vsteps  =  Fidelity v; 

Iplusteps  =  Fidelityl ; 

IpO  *  glpO [ vpump ]  /  glpO [Pumplinecenter ] 

Iplusteps i ze [vpump_ ]  =  - ; 

Iplusteps 

(*Ip0*gIp0  [vpump  ]  /glpO  [  v31  ] 

vlaser  =vhy21  [2,2,1, 87] ; 
vpump  =v31 ; 

o21emsloop =o21ems [vlaser] ; 
o31emsloop =o31ems [vpump] ; 

gapproxf ixtotal211oop =gapproxf ixtotal21 [vlaser] ; 
glp01oop=glp0 [vpump] ; 

Intenfunc=NDSolve  [  {Deqnl  [  z ,  vlaser  ,  vpump  ]  ==Ip  '  [z]  ,  Deqn2  [  z,  vlaser  ,  vpump  ]  ==Iplu  '  [  z]  , 

Deqn3  [  z,  vlaser  ,  vpump  ]  ==Iminu  '  [  z]  ,  Ip  [  znot  ]  ==Ip0*gIp0  [vpump  ]  /  glpO  [  v31  ]  ,  Iplu  [  znot  ]  ==1000 , 
Iplu [znot] ==rl*tl2*Iminu [znot] } , {Ip [z] , Iplu [z] , Iminu [z] } , {z, znot, zfinal} ] *) 
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Percentdone  =  0; 
vpumpcount  =  1 ; 

Do [ j Iplucount  =  1; 

o21emsloop  =  o21ems [vlaser] ; 
a31emsloop  =  a3 1 ems  [vpump] ; 

gapproxf ixtotal211oop  =  gapproxfixtotal21 [vlaser] ; 
glpOloop  =  glpO [vpump] ; 
varray  [  vpumpcount ]  =  vpump ; 

Do[ 

{ 

Intenfunc  =  NDSolve [ {Deqnl [ z ,  vlaser,  vpump]  ==  Ip  '  [z]  , 

Deqn2[z,  vlaser,  vpump]  ==  Iplu'[z],  Deqn3[z,  vlaser,  vpump]  ==  Iminu'[z],  Ip[znot] 
IpO  *  glpO [vpump]  /  glpO [Pumplinecenter ] ,  Iplu[znot]  + Iminu[znot]  ==  Ipluindex, 
Iplu[znot]  ==  rl  *  tl2  *  Iminu  [  znot  ]  }  ,  {Ip,  Iplu,  Iminu},  {z,  znot,  zfinal}]; 
Savelntenfunc [vpumpcount ,  Iplucount ]  =  Intenfunc, 

Iplucount  =  Iplucount  + 1 ; 

_  100.  I  Iplucount -1  \ 

( *  Percentdonez  =Percentdone  + - *  -  , 

vsteps  \  Iplusteps  / 

Print [Percentdone2 , "%  Completed"] *) 

} 

,  {Ipluindex,  0,  IpO  * glpO [vpump ]  /  glpO [Pumplinecenter ] ,  Iplustepsize [vpump ]} j , 

(vpumpcount  -  1 .  ' 

-  , 

vsteps 

Print [Percentdone ,  "%  Completed"]  , 
vpumpcount  =  vpumpcount  + lj 

,  jvpump,  Pumplinecenter  -  vrange  ,  Pumplinecenter  +  vrange  , 


(Pumplinecenter  +  vrange)  -  (Pumplinecenter  - vrange) 

vsteps 


o.% 

Completed 

10.% 

Completed 

20.% 

Completed 

30.% 

Completed 

40.% 

Completed 

50.% 

Completed 

60.% 

Completed 

70.% 

Completed 

80.% 

Completed 
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E.  Routine  for  Selecting  the  Best  Solution  from  NDSolve 

ClearAll [CurrentRoot ,  Rootfindeqn] 

Percentdone  =  0; 

Do[ 

| Rootfindeqn [vpumpcount2 ]  =  IpO, 

Do  [ 

{ (*Print [  Evaluate [Iplu [zfinal] / . Savelntenfunc [vpumpcount2 , Iplucount2 ] ] ==={0 .}],*) 
zarraytest  [vpumpcount2  ,  Iplucount2 ]  =  znot, 

Iparraytest [vpumpcount2  ,  Iplucount2 ]  = 

Evaluate [Ip [zarraytest  [vpumpcount2  ,  Iplucount2 ] ]  /  . 

Savelntenfunc [ vpumpcount2  ,  Iplucount2 ] ]  [ [1] ] , 

Ipnolasingarraytest  [vpumpcount2  ,  Iplucount2 ]  =  Evaluate  [ 

Ip [ zarraytest  [vpumpcount2  ,  Iplucount2 ] ]  / .  Savelntenfunc  [vpumpcount2  ,  1] ]  [ [1] ] , 
Ipluarraytest  [vpumpcount2  ,  Iplucount2 ]  = 

Evaluate [Iplu [zarraytest  [vpumpcount2  ,  Iplucount2 ] ]  /  . 

Savelntenfunc [vpumpcount2  ,  Iplucount2 ] ]  [ [1] ] , 

Iminuarraytest  [vpumpcount2  ,  Iplucount2 ]  = 

Evaluate [Iminu [zarraytest  [vpumpcount2  ,  Iplucount2 ] ]  /  . 

Savelntenfunc [vpumpcount2  ,  Iplucount2 ] ]  [ [1] ] , 

Zarraytest  [vpumpcount2  ,  Iplucount2 ]  = 

Popl [ zarraytest  [vpumpcount2  ,  Iplucount2 ] ,  vlaser ,  varray [vpumpcount2 ] ]  /. 

{ Ip [ zarraytest  [vpumpcount2  ,  Iplucount2 ] ]  ->  Iparraytest  [vpumpcount2  ,  Iplucount2 ] , 

Iplu [ zarraytest  [vpumpcount2  ,  Iplucount2 ] ]  ->  Ipluarraytest [vpumpcount2  ,  Iplucount2 ] , 
Iminu [zarraytest [vpumpcount2  ,  Iplucount2 ] ]  -> 

Iminuarraytest  [vpumpcount2  ,  Iplucount2 ]  }  , 

/Zarraytest  [vpumpcount2 ,  Iplucount2 ]  = 

Pop2 [ zarraytest  [vpumpcount2  ,  Iplucount2 ] ,  vlaser,  varray  [vpumpcount2  ]  ]  /. 

{ Ip [ zarraytest [vpumpcount2  ,  Iplucount2 ] ]  ->  Iparraytest [vpumpcount2 ,  Iplucount2 ] , 

Iplu [ zarraytest [vpumpcount2 ,  Iplucount2 ] ]  ->  Ipluarraytest [vpumpcount2 ,  Iplucount2 ] , 
Iminu [zarraytest [vpumpcount2  ,  Iplucount2 ] ]  -> 

Iminuarraytest [ vpumpcount2  ,  Iplucount2 ]  }  , 

/Zarraytest  [vpumpcount2  ,  Iplucount2 ]  = 

Pop3 [ zarraytest  [vpumpcount2  ,  Iplucount2 ] ,  vlaser,  varray [vpumpcount2 ] ]  /. 

{ Ip [ zarraytest [vpumpcount2  ,  Iplucount2 ] ]  ->  Iparraytest [vpumpcount2  ,  Iplucount2 ] , 

Iplu [ zarraytest [vpumpcount2 ,  Iplucount2 ] ]  ->  Ipluarraytest [vpumpcount2 ,  Iplucount2 ] , 
Iminu [zarraytest [vpumpcount2  ,  Iplucount2 ] ]  -> 

Iminuarraytest [ vpumpcount2  ,  Iplucount2 ] } , 

ytest [vpumpcount2  ,  zarraytest  [vpumpcount2  ,  Iplucount2 ] ]  =  o21ems  [vlaser]  * 
(/Zarraytest  [vpumpcount2 ,  Iplucount2 ]  -  Zarraytest  [vpumpcount2 ,  Iplucount2 ] ) , 

If  [ 
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Evaluate [Iplu [zfinal]  /  .  Savelntenfunc  [vpumpcount2  ,  Iplucount2 ] ] [[1]]  >=  0., 
CurrentRoot [vpumpcount2  ,  Iplucount2 ]  = 

(r2  *  tl A 2  *  Evaluate [ Iplu [ zfinal ]  /  .  Savelntenfunc [vpumpcount2 ,  Iplucount2 ] ]  [ [1] ] )  / 
( (Evaluate [Iminu [zfinal]  /  .  Savelntenfunc [vpumpcount2 ,  Iplucount2 ] ]  [ [1]  ]  )  + 
0.00001)  -  1 

]  , 

If  [ 

ytest  [vpumpcount2  ,  zarraytest [vpumpcount2  ,  Iplucount2 ] ]  -  a  >  0, 

{ If [ Iplucount2  >  1, 

If  [ 

Abs [CurrentRoot [ vpumpcount2 ,  Iplucount2 ] ]  <  Rootfindeqn [vpumpcount2 ] , 

{Rootfindeqn  [vpumpcount2  ]  =  Abs [CurrentRoot  [vpumpcount2  ,  Iplucount2 ] ] , 

Iplubest [ vpumpcount2  ]  =  Iplucount2  } 

] 

]}, 

If [Rootfindeqn [vpumpcount2  ]  ==  IpO ,  {Rootfindeqn  [vpumpcount2  ]  = 

Abs [CurrentRoot [vpumpcount2 ,  1]],  Iplubest [vpumpcount2 ]  =  1}] 

] 

}, 

{Iplucount2,  1,  Iplucount-1,  1} 

]  , 

(vpumpcount2  *  1 . 

— 
vpumpcount  -  1 

Print [Percentdone ,  "%  Completed"] 

(*, Print [Iplubest [vpumpcount2 ] ] *)  ,  Print [Rootfindeqn  [vpumpcount2  ]  ]  j, 

{vpumpcount2  ,  1,  vpumpcount  -1,  1} 

] 

Rootfindavg  = 

1 

-  Sum [ Root findeqn [vpumpcount2 ] ,  {vpumpcount2 ,  1,  vpumpcount  -1,  1 }  ]  ; 

vpumpcount  -  1 

9.09091%  Completed 

1 . 

18.1818%  Completed 

1 . 

27.2727%  Completed 
0.0815172 

36.3636%  Completed 
2 . 15263 

45.4545%  Completed 
3.1977 

54.5455%  Completed 
3.44164 

63.6364%  Completed 
3.1974 
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72.7273%  Completed 
2 . 15352 

81.8182%  Completed 
0.0819762 

90.9091%  Completed 

1 . 

100.%  Completed 
1 . 
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F.  Routine  for  Creating  Arrays  for  Output  Graphs  Based  on  the  Best  Solutions 

zcountfinal  =  Fidelityz ; 

yabovel  -  0; 

Percentdone  =  0; 

Do[ 

{ 

yabovel a  =  0, 

Do[ 

{ 

zfinal  -  znot  zfinal  -  znot 

zarray2  [zcountg]  =  zcountg  *  -  +  znot  -  - , 

zcountfinal  zcountfinal 

zfinal  -  znot  zfinal  -  znot 

zarray  [vpumpcountg  ,  zcountg]  =  zcountg  *  -  +  znot  -  - , 

zcountfinal  zcountfinal 

Iparray [ vpumpcountg  ,  zcountg]  =  Evaluate [ Ip [ zarray [ vpumpcountg  ,  zcountg]]  /. 
Savelntenfunc [vpumpcountg  ,  Iplubest  [vpumpcountg  ]  ]  ]  [  [1] ] , 

Ipnolasingarray [vpumpcountg  ,  zcountg]  = 

Evaluate [Ip [zarray [vpumpcountg ,  zcountg]]  /.  Savelntenfunc [vpumpcountg ,  1]][[1]], 
Ipluarray [ vpumpcountg  ,  zcountg]  =  Evaluate [Iplu [zarray [vpumpcountg  ,  zcountg]]  /. 
Savelntenfunc [vpumpcountg  ,  Iplubest  [vpumpcountg  ]  ]  ]  [  [1] ] , 

Iminuar ray [vpumpcountg  ,  zcountg]  =  Evaluate [ Iminu [ zarray [vpumpcountg ,  zcountg]]  /. 
Savelntenfunc [vpumpcountg  ,  Iplubest  [vpumpcountg  ]  ]  ]  [  [1] ] , 

/larray  [vpumpcountg,  zcountg]  =  Popl [zarray [vpumpcountg ,  zcountg],  vlaser ,  varray [ 
vpumpcountg]]  /.  { Ip [ zarray [vpumpcountg  ,  zcountg]]  ->  Iparray [vpumpcountg  ,  zcountg], 
Iplu [ zarray [ vpumpcountg ,  zcountg]]  ->  Ipluarray  [vpumpcountg  ,  zcountg], 

Iminu [zarray [vpumpcountg  ,  zcountg] ]  ->  Iminuarray [vpumpcountg  ,  zcountg] } , 

/2array  [vpumpcountg,  zcountg]  =  Pop2 [zarray [vpumpcountg ,  zcountg],  vlaser,  varray  [ 
vpumpcountg]]  /.  { Ip [ zarray [vpumpcountg  ,  zcountg]]  ->  Iparray [vpumpcountg  ,  zcountg], 
Iplu [ zarray [ vpumpcountg  ,  zcountg]]  -»  Ipluarray  [vpumpcountg  ,  zcountg], 

Iminu [zarray [vpumpcountg ,  zcountg] ]  -»  Iminuarray [vpumpcountg ,  zcountg] } , 

/Sarray  [vpumpcountg,  zcountg]  =  Pop3 [zarray [vpumpcountg ,  zcountg],  vlaser,  varray  [ 
vpumpcountg]]  /.  { Ip [ zarray  [vpumpcountg  ,  zcountg]]  -*  Iparray [vpumpcountg  ,  zcountg], 
Iplu [ zarray [ vpumpcountg ,  zcountg]]  -»  Ipluarray [vpumpcountg ,  zcountg], 

Iminu [zarray [vpumpcountg  ,  zcountg] ]  -»  Iminuarray  [vpumpcountg  ,  zcountg] } , 
y [ vpumpcountg  ,  zarray [ vpumpcountg  ,  zcountg]]  = 
o21ems  [vlaser  ]  *  (/2array  [vpumpcountg,  zcountg]  -  /iarray  [vpumpcountg,  zcountg]). 

If [yabovela  ==  0,  vypumpcountgmin  =  vpumpcountg] 

If [y [vpumpcountg  ,  zarray [ vpumpcountg  ,  zcountg]]  >  0, 

{yabovel  =  1,  yabovela  =  1,  vjarray  [vpumpcountg,  zcountg]  =  varray [ vpumpcountg ]  , 
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vypumpcountgmax  =  vpumpcountg  }  ]  , 

yarray  [  vpumpcountg  ,  zcountg]  =  y [  vpumpcountg  ,  z array [vpumpcountg  ,  zcountg] ] 

}- 

{  vpumpcountg ,  1 ,  vpumpcount  -  1 ,  1 }  1  , 


Percentdone  =  100  * 

Print [Percentdone , 
If [yabovela  ==  1, 

{ 


' zcountg  *  1 .  ^ 
zcountfinal J 
"%  Completed"], 


IpluTable  =  Table [ {varray [vpumpcountg ] ,  Ipluarray [vpumpcountg  ,  zcountg] } , 
{vpumpcountg,  vjpumpcountgmin  ,  vypumpcountgmax  ,  1}], 

IminuTable  =  Table [ {varray [vpumpcountg ]  ,  Iminuarray [vpumpcountg ,  zcountg] } , 
{vpumpcountg,  vjpumpcountgmin,  vjpumpcountgmax  ,  1}], 

Ipluint  =  Interpolation [IpluTable] , 

Iminuint  =  Interpolation [ IminuTable ] , 

Ipluofz [ zcountg]  = 

NIntegrate [ Ipluint [v]  *glp0[v],  {v,  Min [Table [vjarray  [vpumpcountg  ,  zcountg], 
{vpumpcountg,  vjpumpcountgmin  ,  vjpumpcountgmax  ,  1}]],  Max [Table [vjarray  [ 
vpumpcountg,  zcountg],  {vpumpcountg,  vjpumpcountgmin,  vypumpcountgmax,  1}]]}], 
Iminuofz [zcountg]  =  NIntegrate [Iminuint [v]  *  glpO [v] , 

{v,  Min [Table [vjarray  [vpumpcountg  ,  zcountg],  {vpumpcountg,  vjpumpcountgmin  , 
vjpumpcountgmax,  1}]],  Max [ Table [  vjarray  [  vpumpcountg  ,  zcountg], 

{vpumpcountg,  vjpumpcountgmin  ,  vypumpcountgmax  ,  1}] ]  }]  , 
yTable  =  Table  [  {varray [vpumpcountg  ]  ,  y [  vpumpcountg  ,  zarray [vpumpcountg  ,  zcountg]]}, 
{vpumpcountg,  vjpumpcountgmin  ,  vypumpcountgmax,  1}], 
yint  =  Interpolation  [yTable] ,  yofz [zcountg]  =  NIntegrate [yint [v]  *  glpO [v] , 

{v,  Min [Table  [vyarray  [vpumpcountg  ,  zcountg],  {vpumpcountg,  vjpumpcountgmin  , 
vypumpcountgmax,  1}]],  Max [ Table [ vyarray  [  vpumpcountg  ,  zcountg], 

{vpumpcountg,  vjpumpcountgmin  ,  vypumpcountgmax,  1}]]}] 

} 


] 

}■ 

{zcountg,  1,  zcountfinal,  1 } J  (*If  you  really  need  speed  here  comment  out  the  sections 
on  gain  and  their  plotting  routines.  But  you  will  lose  alot  of  information.*) 


10.%  Completed 
20.%  Completed 
30.%  Completed 
40.%  Completed 
50.%  Completed 
60.%  Completed 
70.%  Completed 
80.%  Completed 
90.%  Completed 
100.%  Completed 
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G.  Tabular  and  Plotting  Routines  for  Best  Solutions 

Ipplot  = 

Table [ { varray [vpumpcountg  ]  ,  zarray [vpumpcountg  ,  zcountg] ,  Iparray [ vpumpcountg  ,  zcountg] } , 
{ vpumpcountg  ,  1,  vpumpcount  -1,  1},  {zcountg,  1,  zcountfinal ,  1}]; 

Ipnolasingplot  =  Table [ {varray [vpumpcountg ]  ,  zarray [ vpumpcountg  ,  zcountg], 

Ipnolasingarray [vpumpcountg  ,  zcountg] } , 

{vpumpcountg,  1,  vpumpcount  -1,  1},  {zcountg,  1,  zcountfinal,  1}]; 

Ipluplot  =  Table [ {varray [vpumpcountg  ]  ,  zarray [ vpumpcountg  ,  zcountg], 

Ipluarray [ vpumpcountg  ,  zcountg]}, 

{vpumpcountg,  1,  vpumpcount  -1,  1},  {zcountg,  1,  zcountfinal,  1}]; 

Iminuplot  =  Table [ {varray [vpumpcountg ]  ,  zarray [vpumpcountg  ,  zcountg] , 

Iminuarray [vpumpcountg  ,  zcountg] } , 

{vpumpcountg,  1,  vpumpcount  -1,  1},  {zcountg,  1,  zcountfinal,  1}]; 
yplot  =  Table [ {varray [vpumpcountg ]  ,  zarray [vpumpcountg  ,  zcountg],  yarray [  vpumpcountg , 
zcountg]},  {vpumpcountg,  1,  vpumpcount  -1,  1},  {zcountg,  1,  zcountfinal,  1}]; 
aplot  =  Table [{ varray  [  vpumpcountg  ]  ,  zarray [vpumpcountg  ,  zcountg],  a}, 

{vpumpcountg,  1,  vpumpcount  -1,  1},  {zcountg,  1,  zcountfinal,  1}]; 

Ilaseoutplot  =  Table [{ varray [ vpumpcountg ]  ,  (1  -  r2)  * Ipluarray [vpumpcountg ,  zcountfinal]}, 

{ vpumpcountg  ,  1 ,  vpumpcount  -  1 ,  1 }  ]  ; 

/Iplot  =  Table [{ varray  [  vpumpcountg  ]  ,  zarray [vpumpcountg  ,  zcountg],  flarray  [vpumpcountg, 
zcountg]},  {vpumpcountg,  1,  vpumpcount  -1,  1},  {zcountg,  1,  zcountfinal,  1}]; 

;2plot  =  Table [{ varray [ vpumpcountg  ]  ,  zarray [vpumpcountg  ,  zcountg],  /2array  [vpumpcountg, 
zcountg]},  {vpumpcountg,  1,  vpumpcount  -1,  1},  {zcountg,  1,  zcountfinal,  1}]; 

/Splot  =  Table [{ varray  [  vpumpcountg  ]  ,  zarray [vpumpcountg  ,  zcountg],  /Sarray  [vpumpcountg, 
zcountg]},  {vpumpcountg,  1,  vpumpcount  -1,  1},  {zcountg,  1,  zcountfinal,  1}]; 

Inv31plot  =  Table [{ varray [ vpumpcountg  ]  ,  zarray [vpumpcountg  ,  zcountg], 

/Sarray  [vpumpcountg,  zcountg]  -  /iarray  [vpumpcountg,  zcountg]}, 

{vpumpcountg,  1,  vpumpcount  -1,  1},  {zcountg,  1,  zcountfinal,  1}]; 

Inv21plot  =  Table [{ varray [ vpumpcountg  ]  ,  zarray  [vpumpcountg  ,  zcountg], 

J2array  [vpumpcountg,  zcountg]  - /larray  [vpumpcountg,  zcountg]}, 

{vpumpcountg,  1,  vpumpcount  -1,  1},  {zcountg,  1,  zcountfinal,  1}]; 
yofzplot  =  Table [{ zarray2 [ zcountg] ,  yofz [zcountg] } ,  {zcountg,  1,  zcountfinal,  1}]; 
Ipluofzplot  =  Table [{ zarray2 [ zcountg] ,  Ipluofz [zcountg] } ,  {zcountg,  1,  zcountfinal,  1}]; 
Iminuofzplot  =  Table [{ zarray2 [ zcountg] ,  Iminuofz [zcountg] } ,  {zcountg,  1,  zcountfinal,  1}]; 

ListPlot3D [ Ipplot,  DataRange  ->  All,  Mesh  ->  25, 

InterpolationOrder  ->  All,  ColorFunction  -»  "SouthwestColors " , 

PlotLabel  ->  Style [ "Pump  Intensity  vs.  Freq .  vs.  Pos . " ,  24], 

AxesLabel  -»  \  Pump  Freq  [Hz]  ,  Pos  [m]  ,  Ip  —  1,  ( *LabelStyle-»Directive  [Large  ],  *) 
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PlotRange  -»  {{Pumplinecenter  - vrange ,  Pumplinecenter  +  vrange  }  ,  {0,  celllength},  {0,  Ip0}}j 
ListPlot3D  [ipnolasingplot ,  DataRange  -»  All,  Mesh  -*  25, 

InterpolationOrder  ->A11,  ColorFunction  ->  "SouthwestColors  " , 

PlotLabel  ->  Style  [  "Pump  Intensity  (no  lasing)  vs.  Freq.  vs.  Pos . " ,  24], 

AxesLabel  -»  jpump  Freq [Hz] ,  Pos [m] ,  Ip  [ — jj,  (*LabelStyle-»Directive [Large]  ,  *) 

L  m2 

PlotRange  -»  {{Pumplinecenter  -  vrange,  Pumplinecenter  +  vrange  }  ,  {0,  celllength},  {0,  Ip0}}j 
ListPlot3D  [ipluplot,  DataRange  ->  All,  Mesh  ->  None, 

InterpolationOrder  -»  All,  ColorFunction  ->  "SouthwestColors", 

PlotLabel  -»  Style ["Plus  Wave  Intensity  vs.  Freq.  vs.  Pos.",  24], 

W 

AxesLabel  ->  jpump  Freq  [Hz]  ,  Pos  [m]  ,  I+ [ — j  J  (*,  LabelStyle ^Directive  [Large]  *)  j 

ListPlot3D [iminuplot,  DataRange  ->  All,  Mesh  ->  None, 

InterpolationOrder  ->  All,  ColorFunction  -»  "SouthwestColors", 

PlotLabel  ->  Style [ "Minus  Wave  Intensity  vs.  Freq.  vs.  Pos.",  24], 

W 

AxesLabel  jpump  Freq [Hz] ,  Pos [m] ,  I_  [ — j J (*, LabelStyle ^Directive [Large] *) j 

L  m2 

ListPlot3D  J/Iplot  ,  DataRange  ->  All,  Mesh ->  None,  InterpolationOrder  ->A11, 

ColorFunction  -»  "SouthwestColors",  PlotLabel  ->Style["Ni  vs.  Freq.  vs.  Pos.",  24], 

AxesLabel  -»  \  Pump  Freq [Hz] ,  Pos [m] ,  Cone  —  [(*, LabelStyle ^Directive [Large] *) 

L  L  m3  J  J  J 

ListPlot3D  J/2plot  ,  DataRange  ->  All ,  Mesh ->  None,  InterpolationOrder  ->A11, 

ColorFunction  -»  "SouthwestColors",  PlotLabel  ->Style["N2  vs.  Freq.  vs.  Pos.",  24], 

AxesLabel  jpump  Freq  [Hz]  ,  Pos [m] ,  Cone [ — j J  (*, LabelStyle ^Directive [Large] *) j 

L  m3 

ListPlot3D  J/Splot  ,  DataRange  ->  All ,  Mesh  ->  None,  InterpolationOrder  ->A11, 

ColorFunction  -»  "SouthwestColors",  PlotLabel  ->Style["N3  vs.  Freq.  vs.  Pos.",  24], 

AxesLabel  -»  jpump  Freq [Hz] ,  Pos [m] ,  Cone [ — j J  (*, LabelStyle ^Directive [Large] *) j 

ListPlot3D  [lnv31plot,  DataRange  ->  All,  Mesh  ->  None,  InterpolationOrder  ->  All, 

ColorFunction  -»  "SouthwestColors",  PlotLabel  ->  Style [ "N3-N1  vs.  Freq.  vs.  Pos.",  24], 

AxesLabel  ->  jpump  Freq  [Hz]  ,  Pos [m] ,  Cone [ — j J  (*, LabelStyle ^Directive [Large] *) j 

L  m3 

ListPlot3D [lnv21plot ,  DataRange  ->  All,  Mesh  ->  None,  InterpolationOrder  -»  All, 

ColorFunction  -»  "SouthwestColors",  PlotLabel  ->  Style [ "N2-N!  vs.  Freq.  vs.  Pos.",  24], 
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AxesLabel  ->  jpump  Freq  [Hz]  ,  Pos  [m]  ,  Cone  [ — JJ(*, 

L  m3 

PlotRange->{  Automatic , Automatic ,{-1017,1017}}*)  (* , LabelStyle-»Directive [Large] *) j 

ListPlot3D  [  {yplot ,  aplot  }  ,  DataRange  ->  All,  Mesh ->  None,  ColorFunction  -»  "SouthwestColors  "  , 
InterpolationOrder  -»A11,  PlotLabel  -»  Style  [  "Gain  Coeffcient  vs.  Freq.  vs.  Pos.",  24], 
AxesLabel  ->  {Pump  Freq  [Hz]  ,  Pos  [m]  ,  y}  , 

PlotRange  ->  {  {Pumplinecenter  -  vrange  ,  Pumplinecenter  +  vrange  }  ,  {0,  celllength},  {-2a,  2  a}} 
(*, LabelStyle-»Directive [Large] *) ] 

ListLinePlot [ {yof zplot } ,  PlotLabel  -> Style [ "y  (Integrated)  vs.  Pos",  24], 

AxesLabel  ->  {Pos [m] ,  y} (*, LabelStyle-»Directive [Large] *) ] 

ListLinePlot [ { Ipluof zplot ,  Iminuof zplot } , 

PlotLabel  ->Style["I  +  and  I_ (Integrated  over  vp)  vs.  Pos",  24], 

AxesLabel  ->  jpos [m] ,  I [ — J J (*, LabelStyle-»Directive [Large] *) ] 

ListLinePlot [ilaseoutplot ,  PlotLabel  ->  Style [ "Laser  Output  Intensity",  24], 

W 

AxesLabel  ->  jpump  Freq  [Hz]  ,  Iout  [ — j  J  (*,  LabelStyle ^Directive  [Large ]  *)  j 
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H.  Laser  Output  Power  and  Other  System  Outputs 

Laseroutfunc  =  Interpolation [ Ilaseoutplot ] ; 

Plot [Laseroutfunc [ v]  *  glpO [ v] ,  {  v,  Pumplinecenter  -  vrange ,  Pumplinecenter  +  vrange  }  ]  ; 
LasPower  =  NIntegrate [Laseroutfunc [v]  *  glpO [v] , 

{v,  Pumplinecenter  -  vrange  ,  Pumplinecenter  +  vrange } ] ; 
yofzint  =  Interpolation [yofzplot ]  ; 

1 

ytotal  =  - NIntegrate  [yofzint  [z]  ,  {z,  znot,  zfinal}]; 

zfinal  -  znot 

OutputEf f iciency  =  LasPower  /  IpO ; 

(♦Interpolation  [  yplot ] 

If [yabovel =  =  1 , ygain  =NIntegrate [ 

NIntegrate [ Interpolation [  yplot ] *glp0 [ v] , { v, Min [ v-jarray  [ vpumpcountg , zcountg] ] , 

Max  [vnarray  [vpumpcountg  ,  zcountg]  ]  }  ]  ,  {z,  znot,  zfinal}  ]  ]  ;  *) 

Print [style [ "1 .  The  Output  Laser  Intensity  of  DPAL  System  Tested  is  predicted  to  be  ", 
30,  Bold],  Style [LasPower ,  30,  Bold,  Red],  Style [ "W/m2  "  ,  30,  Bold,  Red]] 

Print [Style [ "2 .  The  Output  Lasing  Frequency  of  DPAL  System  Tested  is  predicted  to  be  ", 
30,  Bold],  Style [vlaser ,  30,  Bold,  Red],  Style["Hz",  30,  Bold,  Red]] 

Print [Style [ "3 .  The  Output  Efficiency  of  the  Tested  DPAL  System  is  predicted  to  be  ", 
30,  Bold],  Style [OutputEf ficiency  *  100 ,  30,  Bold,  Red],  Style["V,  30,  Bold,  Red]] 

Print [style [ "4 .  Ip0/Isat  for  this  system  is  ",  30,  Bold], 

Style [ IpO  /  Isat [Pumplinecenter ] ,  30,  Bold,  Red]] 

Print [Style [ "5 .  y(total)/a  or  Gain  Divded  by  Loss  for  this  system  is  ",  30,  Bold], 

Style [ytotal  /a,  30,  Bold,  Red]] (*This  is  a  measure  of  the  accuracy  of  the  simulation*) 
Print [Style [ "6 .  The  degree  to  which  the  boudary  condition  at  z  final  was  met  was:  ", 

30,  Bold],  Style [Rootfindavg ,  30,  Bold,  Red]] (*, 

"Note  The  smaller  this  number  is  the  better .  Above  1  is  bad  and  implies 
a  solution  set  should  be  rerun  with  a  longer /shorter  cavity.  "*) 

(♦This  is  a  measure  of  the  accuracy  of  the  simulation*) 
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1 .  The  Output  Laser  Intensity 
of  DPAL  System  Tested  is 

predicted  to  be  358  130.W/m2 

2 .  The  Output  Lasing  Frequency 
of  DPAL  System  Tested  is 

predicted  to  be  3 . 77109  x  1014Hz 

3 .  The  Output  Efficiency 
of  the  Tested  DPAL  System 
is  predicted  to  be  0.71626% 

4.  Tp0/Tsat  for  this  system  is 
103.782 

5.  y (total) /a  or  Gain  Divded  by 
Loss  for  this  system  is  154.126 

6 .  The  degree  to  which 

the  boudary  condition  at  z 
final  was  met  was:  1.66422 
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